Imaging a subsurface geological model at a past intermediate restoration time

ABSTRACT

A system and method is provided for restoring a 3D tomographic model of the Earth&#39;s subsurface geology from the present-day to a past restoration time. Whereas at the present time all faults represent active discontinuities, at a past restoration time some faults have not yet formed. Accordingly, the restored model divides the fault network into τ-active faults (discontinuous surfaces for faults that intersect the layer deposited at the past restoration time) and τ-inactive faults (continuous surfaces for faults that do not intersect the layer deposited at the past restoration time). A new 3D restoration transformation is also provided that uses linear geological constraints to process the restoration model in less time and generate more accurate geological images.

CROSS REFERENCE TO RELATED APPLICATIONS

This application is a continuation-in-part of U.S. Ser. No. 16/681,061,filed on Nov. 12, 2019, which in turn is a continuation of U.S. Ser. No.16/244,544, filed on Jan. 10, 2019, both of which are incorporatedherein by reference in their entirety.

FIELD OF THE INVENTION

Embodiments of the invention relate to the field of geologicaltomography for generating an image of the interior subsurface of theEarth based on geological data collected by transmitting a series ofincident waves and receiving reflections of those waves acrossdiscontinuities in the subsurface. The incident and reflected waves arereconstituted by a 3D model to generate an image of the reflectingsurfaces interior to the Earth. Accordingly, geological tomographyallows geophysicists to “see inside” the Earth.

Embodiments of the invention further relate to geological restoration inwhich the tomographic images of the present day geology are transformedinto images of the past geology, as it was configured at an intermediaterestoration time in the past τ before the present day and after thestart of deposition of the oldest subsurface layer being imaged. Newtechniques are proposed herein to improve both the accuracy andcomputational speed of generating those images of the past restoredgeology. Improved images may aid geoscientists exploring the subsurfacegeology for applications such as predicting tectonic motion orearthquakes, or by engineers in the mining or oil and gas industries.

BACKGROUND OF THE INVENTION

The accuracy of a geological model of the present day configuration ofthe subsurface of the Earth may be improved by “restoring” the model toa past intermediate time τ and checking model consistency at that timein the past. However, restoring geological models is a complex task andcurrent methods are typically inefficient, requiring extensiveprocessing resources and time, as well as inaccurate, relying onover-simplifications that induce errors to moderate the complexity ofthe task.

There is a longstanding need in the art to efficiently and accuratelyrestore geological models from their present day geology to their pastgeology at restored past time τ.

SUMMARY OF EMBODIMENTS OF THE INVENTION

Some embodiments of the invention are directed to modeling restoredgeological models with τ-active and τ-inactive faults. In an embodimentof the invention, a system and method is provided for restoring a 3Dmodel of the subsurface geology of the Earth from a present day geometrymeasured at a present time to a predicted past geometry at a pastrestoration time. The 3D model of the present day measured geometrycomprising a network of faults may be received, wherein a fault is adiscontinuity that divides fault blocks that slide in oppositedirections tangential to the surface of the fault as time approaches amodeled time. A past restoration time τ may be selected that is prior tothe present time and after a time when an oldest horizon surface in the3D model was originally deposited. The network of faults may be dividedinto a subset of τ-active faults and a subset of τ-inactive faults,wherein a τ-active fault is a fault that is active at the pastrestoration time τ and a τ-inactive fault is a fault that is inactive atthe past restoration time τ. A fault may be determined to be τ-activewhen the fault intersects a horizon H_(τ) that was originally depositedat the past restoration time τ and a fault may be determined to beτ-inactive when the fault does not intersect the horizon H_(τ) that wasoriginally deposited at the past restoration time τ. The 3D model may berestored from the present day measured geometry to the predicted pastgeometry at the past restoration time τ by modeling each τ-active andτ-inactive fault differently. Each τ-active fault may be modeled to joinend points of a horizon H_(τ) separated on opposite sides of the faultin the present day model to merge into the same position in the restoredmodel by sliding the end points towards each other in a directiontangential to the surface of the τ-active fault. Each τ-inactive faultmay be modeled to keep collocated points on opposite sides of the faulttogether.

Some embodiments of the invention are directed to modeling restoredgeological models with new restoration coordinates u_(τ), v_(τ), t_(τ).In an embodiment of the invention, a system and method is provided forrestoring a 3D model of the subsurface geology of the Earth from apresent day measured geometry to a predicted past geometry at arestoration time in the past τ. The 3D model of the present day geometryof the subsurface may be received, including one or more foldedgeological horizon surfaces. A value may be selected of a restorationtime in the past τ before the present day and after a time an oldesthorizon surface in the 3D model of the subsurface was deposited. The 3Dmodel may be restored from the present day measured geometry to thepredicted past geometry at the restoration time in the past τ using a 3Dtransformation. The vertical component of the 3D transformation mayrestore the geometry to the vertical coordinate t_(τ) such that: pointsalong a horizon surface H_(τ) modeling sediment that was deposited atthe selected restoration time in the past τ have a substantiallyconstant value for the restored vertical coordinate t_(τ); and at anylocation in the 3D model, the restored vertical coordinate t_(τ) isequal to a sum of a first approximation t′_(τ) of the verticalcoordinate and an error correction term ε_(τ), wherein the errorcorrection term ε_(τ) is computed by solving a linear relationship inwhich a variation in the sum of the first approximation t′_(τ) of thevertical coordinate and the error correction term ε_(τ) between any twopoints separated by an infinitesimal difference in the direction ofmaximal variation of the sum is approximately equal to the distancebetween the points in the direction of maximal variation; and displayingan image of the restored 3D model of the subsurface geology of the Earthsuch that each point in the 3D model is positioned at the restoredvertical coordinate t_(τ) as it was configured at the restoration timein the past τ.

Some embodiments of the invention are directed to modeling restoredgeological models taking compaction into account at an intermediaterestoration time in the past τ. In an embodiment of the invention, asystem and method is provided for decompacting a 3D model of thesubsurface geology of the Earth at an intermediate restoration time inthe past τ. Some embodiments may receive a 3D model of present-daygeometry of the subsurface geology and a measure of present-day porosityexperimentally measured within the subsurface geology of the Earth. Avalue of a restoration time in the past τ may be selected before thepresent day and after a time an oldest horizon surface in the 3D modelof the subsurface was deposited. The 3D model from the present daymeasured geometry may be restored to the predicted past geometry at therestoration time in the past τ using a 3D transformation. The verticaldimension of the restored 3D model may be decompacted to elongatevertical lengths of geological layers below a horizon layer deposited atthe restoration time in the past τ. The vertical lengths may beelongated based on a relationship between a depositional porosity of thegeological layers at the time sediment in those layers was deposited,restoration porosity of the geological layers at the restoration time inthe past τ, and the present-day porosity of the geological layersexperimentally measured in the present-day.

BRIEF DESCRIPTION OF THE DRAWINGS

The principles and operation of the system, apparatus, and methodaccording to embodiments of the present invention may be betterunderstood with reference to the drawings, and the followingdescription, it being understood that these drawings are given forillustrative purposes only and are not meant to be limiting.

FIG. 1 schematically illustrates an exploded view of a 3D present daygeological model of the subsurface of the Earth according to anembodiment of the invention. The 3D geological model may comprise afaulted 3D grid Γ 100. Cell edges 106 of the grid are constrained tonever cross faults 105. During restoration, twin faces F⁻ 104 and F⁺ 103on opposite sides of a fault F 105 may slide along one another only in adirection tangential to the surface of the fault F 105. Points (r_(F) ⁺,τ_(F) ⁻) (101, 102) are twin-points. “Twin” points or faces may refer topoints or faces that were collocated at the time of their deposition,but which may have separated at a later time.

FIG. 2 schematically illustrates a vertical cross section of a 3Dsubsurface model in a volume deformation: the u_(τ)v_(τ)t_(τ)-transform201 restores the volume G_(τ) 202 in a present day space G 220 to a pastvolume G _(τ) 203 in a past restored space 219 as the subsurface wasconfigured at restoration time τ according to an embodiment of theinvention.

FIG. 3 schematically illustrates a vertical cross section of the presentday domain G_(τ) 202 according to an embodiment of the invention.Horizons 210 are shown in bold lines and level sets 208 of the verticalrestoration coordinate t_(τ)(r) are shown in dashed lines. τ-activefaults 105 (depicted as bold black lines) cut the horizon H_(τ) 210deposited at restoration time_(τ) whilst τ-inactive faults 300 (depictedas bold gray lines) do not cut horizon H_(τ) 210. Level sets of thevertical restoration coordinate t_(τ)(r) (dashed lines 208) arecontinuous across τ-inactive faults 300 and only separated acrossτ-active faults 105 because the sedimentary layers of those level setswere not divided when the faults were inactive (at a time before thefaults formed).

FIG. 4 schematically illustrates an example 3D geological model used totest the accuracy of the vertical restoration coordinate t_(τ)(r)according to an embodiment of the invention. This example 3D geologicalmodel has a “ramp” geological structure with a restored horizon H_(τ)210 as the central, sigmoid surface. This seemingly “simple” testcomprises highly complex calculations because there are large variationsin layer thickness which make numerical computation of the geologicaltime of deposition t(r) more challenging.

FIG. 5 shows a comparison of multiple histograms 501, 502 and 503 of∥grad t_(τ)(r)∥ each using a different method to compute t_(τ)(r) in thetest case domain G_(τ) of FIG. 4. Depending on the method used tocompute t_(τ)(r), the resulting magnitude of grad t_(τ)(r) may severelydeviate from the ideal value of “1,” as required by equation (10).Histogram 503, which is generated according to an inventive embodiment,is the closest approximation of (10), yielding the most accuraterestoration model.

FIG. 6 schematically illustrates a vertical cross section of a 3Dsubsurface model in which a pair of τ-twin-points ({tilde over (r)}_(F)⁺, r_(F) ⁻)_(τ) (601, 602) are deduced from a pair of present daytwin-points (r_(F) ⁺, r_(F) ⁻) (101, 102), using a depositional (e.g.,GeoChron) model as input, according to an embodiment of the invention.

FIG. 7 schematically illustrates a direct uvt-transform 700 and inverseuvt-transform 701 that transform a 3D model between a present daygeological space G 220 and a depositional geological space G 719according to an embodiment of the invention.

FIG. 8 schematically illustrates a u_(τ)v_(τ)t_(τ)-transform 201 of avertical cross-section of a 3D subsurface model showing the geologicalimpact between fault 300 being erroneously considered a τ-active fault(top image of FIG. 8) vs. correctly considered a τ-inactive fault withrespect to restored horizon H_(τ) 210 (bottom image of FIG. 8). In thetop image of FIG. 8, when fault 300 is erroneously considered a τ-activefault, in order to preserve geological volume, present day fault block800 is transformed to a restored fault block 801 that intersectsτ-active fault 805, which contradicts geological rules. In contrast, inthe bottom image of FIG. 8, when fault 300 is correctly considered aτ-inactive fault, a volume-preserving transformation maps the presentday fault block 800 correctly to a restored fault block that stayswithin (and does not cross) τ-active fault 805, according to geologicalrules.

FIG. 9 shows a comparison of multiple histograms 901, 902 and 903 ofΔV/V each using a different method to compute t_(τ)(r) in the test casedomain G_(τ) of FIG. 4. Depending on the method used to computet_(τ)(r), the resulting magnitude of ΔV/V may severely deviate from thetarget value of “0”. Histogram 903, which is generated according to anembodiment of the invention, has the least amount of volume variation,yielding the most accurate restoration model.

FIG. 10 schematically illustrates an example of a 3Du_(τ)v_(τ)t_(τ)-transform 201 of a horizon H_(τ) 210 (depicted as awhite layer) from a present day domain G 202 to a restored domain G _(τ)203 at restoration time τ according to an embodiment of the invention.

FIG. 11 schematically illustrates an example of a sequence ofrestorations at geological times {τ₁<τ₂< . . . <τ_(n)} of vertical crosssections of the 3D model of FIG. 10 according to an embodiment of theinvention.

FIG. 12 schematically illustrates an example of fault striae induced onfaults by paleo-geographic coordinates u(r) and v(r) of a depositionalmodel provided as input to the restoration according to an embodiment ofthe invention.

FIG. 13 schematically illustrates an example display of a sequentialchronological transformation of a vertical cross section of a 3Dsubsurface model according to an embodiment of the invention.Chronological transformation may correspond to a plurality ofrestoration times τi ordered in a sequence either “forward modeling”according to a forward passage of time (e.g., with sequentiallyascending or later values of time) or “restoration” according to areverse passage of time (e.g., with sequentially descending or earliervalues of time). In the example shown in FIG. 13, the sequence ofchronological transformation progresses according to a forward passageof time, from the start of deposition at earliest geological time τ₁, toone or more subsequent intermediate restorations times τ₂<τ₃, to alatest present day geological time τ₄ (although the sequence mayprogress according to the reverse passage of time). Further, any othernumber or orders of restorations or times may be used).

FIG. 14 schematically illustrates a geological tomography technique inwhich a series of incident and reflected waves are propagated through asubsurface region of the Earth to image the subsurface according to anembodiment of the invention.

FIG. 15 schematically illustrates a system for restoring a present daygeological model to an intermediate restoration time τ, according to anembodiment of the invention.

FIG. 16 is a flowchart of a method to restore a geological model withτ-active and τ-inactive faults, according to an embodiment of theinvention.

FIG. 17 is a flowchart of a method to restore a geological model withimproved accuracy using a new thickness-preserving constraint, accordingto an embodiment of the invention.

FIG. 18 schematically illustrates an example 3D geological volume of amodel that is compacted (right image) and that is decompacted (leftimage) at an intermediate restoration time in the past τ, according toan embodiment of the invention.

FIG. 19 is a flowchart of a method for decompacting a 3D model of thesubsurface geology of the Earth at an intermediate restoration time inthe past τ, according to an embodiment of the invention.

For simplicity and clarity of illustration, elements shown in thedrawings have not necessarily been drawn to scale. For example, thedimensions of some of the elements may be exaggerated relative to otherelements for clarity. Further, where considered appropriate, referencenumerals may be repeated among the drawings to indicate corresponding oranalogous elements throughout the serial views.

DETAILED DESCRIPTION OF EMBODIMENTS OF THE INVENTION

Embodiments of the invention improve conventional restoration techniquesfor imaging restored geological models as follows:

“τ-Active” Faults Vs. “τ-Inactive” Faults:

In conventional restoration models, all faults are active (asdiscontinuous surfaces) at all times. However, in reality, certainfaults have not yet formed or activated at various intermediaterestoration times τ. Accordingly, conventional restoration modelsgenerate false or “phantom” faults that erroneously divide geology thathas not yet fractured, leading to geological inaccuracies in subsurfaceimages.

Embodiments of the invention solve this problem by selectivelyactivating and deactivating individual fault surfaces to bediscontinuous or continuous, respectively, depending on the specificrestoration geological-time τ. For each intermediate restoration time inthe past τ, embodiments of the invention split faults into twocomplementary subsets of “τ-active” faults and “τ-inactive” faults.τ-active faults are activated at restoration time τ (e.g., adiscontinuous fault surface along which fault blocks slidetangentially), whereas τ-inactive faults are deactivated at restorationtime τ (e.g., a continuous surface that does not behave as a fault).

As faults form and evolve over time, they behave differently atdifferent geological times in the past. For example, a fault that formsat an intermediate geological-time τ, where τ₁<τ<τ₂, is τ-active in arestored model at later time τ₂ (after the fault has formed), butτ-inactive in a restored model at earlier time τ₁ (before the fault hasformed). This fault classification allows faults to be modelleddifferently at each restoration time τ in a geologically consistent way,thereby preventing unrealistic deformations from being generated in theneighborhood of these faults.

FIG. 8 shows the problem of a fault 300 being erroneously consideredactive at a restoration time before it formed (top image of FIG. 8) andthe solution of modeling the fault as a τ-inactive fault to correctlydeactivate the fault at restoration time τ according to embodiments ofthe invention (bottom image of FIG. 8). In the top image of FIG. 8, whena horizon H_(τ) 210 is restored using u_(τ)v_(τ)t_(τ)-transform 201,fault block 800 (shaded region in the top-left image of FIG. 8) isbounded by an active fault 105 and an inactive fault 300. If howeverfault 300 is erroneously considered as an active fault then, afterapplying restoration 201:

-   -   erroneous τ-twin points (803, 823) are transformed into a pair        of collocated points 813=833,    -   τ-twin points (804, 824) are transformed into collocated points        814=834.    -   It is clear that, if d(a, b) denotes the distance between any        arbitrary pair of points (a, b), then:

d(803,824)=d(833,835)≠d(833,834)   (2)

-   -   This observation shows that erroneously considering fault 300 as        a τ-active fault inevitably generates unrealistic deformations.

This problem is solved according to embodiments of the invention, e.g.,as shown in the bottom image of FIG. 8. In this image, fault 105 ismodeled as a τ-active fault (activating the fault), but fault 300 ismodeled as a τ-inactive fault (deactivating the fault). Accordingly,when u_(τ)v_(τ)t_(τ)-transform 201 is applied to fault block 800 (shadedregion in the bottom-left image of FIG. 8), restored fault block 801(shaded region in the bottom-right image of FIG. 8) is no longer boundedby an active fault (τ-inactive fault 300 is inactive at restored timeτ). Accordingly, the restored fault block 801 preserves volume and stayswithin (and does not cross) τ-active fault 805 (because the deactivatedboundary transformed from τ-inactive fault 300 may shift to accommodatea shift in the restored τ-active fault 805).

Contrary to conventional methods, the use of τ-active and τ-inactivefaults produces more accurate results, e.g., even if there is nocontinuous path between (no way to continuously connect) a given faultblock (e.g., 800) and the horizon H_(τ) (e.g., 210) deposited atgeological time τ, which typically requires additional processing thatmay induce errors. By selectively activating and inactivating faults atthe various restoration times according to when they form, embodimentsof the invention eliminate erroneous phantom faults and more accuratelyrepresent the faulted geology.

Reference is made to FIG. 16, which is a flowchart of a method torestore a geological model using τ-active and τ-inactive faults,according to an embodiment of the invention.

In operation 1610, a processor may receive a 3D model of the present daymeasured geometry comprising a network of faults (e.g., present daymodel 202). The present day model may be measured tomographically byscanning the Earth's subsurface e.g., as described in reference to FIGS.14 and 15. In the present day geology, all faults in the model havealready formed and so, represent active discontinuities that dividefault blocks which slide in opposite directions tangential to thesurface of the fault as time approaches a modeled time.

In operation 1620, a processor may select or receive a past restorationtime τ that is “intermediate” or prior to the present time and after thestart of the subsurface's deposition (the time period when an oldesthorizon surface in the 3D model was originally deposited).

In operation 1630, a processor may divide the network of faults into asubset of τ-active faults and a subset of τ-inactive faults. τ-activefaults may be faults that are active at the past restoration time τ andτ-inactive faults are faults that are inactive at the past restorationtime τ. A fault is determined to be τ-active when the fault intersects ahorizon H_(τ) that was originally deposited at the past restoration timeτ (e.g., see τ-active faults 105 of FIG. 3) and a fault is determined tobe τ-inactive when the fault does not intersect the horizon H_(τ) thatwas originally deposited at the past restoration time τ (e.g., seeτ-inactive faults 300 of FIG. 3). Because different faults activate tofracture the subsurface at different geological times, the processor maydivide the fault network differently at different geological times.Thus, a fault may be τ-active at a first restoration time τ′ (e.g., atime period during which the fault has formed) and τ-inactive at asecond restoration time τ″ (e.g., a time period different than thatduring which the fault has formed). In one embodiment, iso-valuesurfaces (e.g., 208 of FIG. 3) of each restoration coordinate (e.g.,u_(τ), v_(τ), and t_(τ)) are continuous across τ-inactive faults (e.g.,300 of FIG. 3) and discontinuous across τ-active faults (e.g., 105 ofFIG. 3).

In operation 1640, a processor may restore the 3D model from the presentday measured geometry to the predicted past geometry at the pastrestoration time τ. During restoration, the processor may flatten ahorizon H_(τ) (e.g., 210 of FIG. 4) that was originally deposited attime τ to a substantially planar surface of approximately constantdepth. For horizons older (e.g., deposited deeper in the subsurface)than horizon H_(τ), the processor may restore the horizons to non-planarsurfaces, e.g., when the thickness of the layers is not constant.Because the region of the subsurface deposited after the restorationtime τ (e.g., deposited shallower in the subsurface) did not yet existat the time of the restored model, restoring the 3D model to a pastrestoration time τ may eliminate (e.g., removing or not displaying) allrelatively shallower horizon surfaces that were originally depositedafter the past restoration time τ. During restoration, the processor maytreat τ-active and τ-inactive faults differently in operations 1650 and1660, respectively.

In operation 1650, for each τ-active fault, a processor may model theτ-active fault as an active discontinuous fault surface and restore thehorizon surface by removing or omitting the fault surface at the time ofrestoration. The processor may eliminate the τ-active fault duringrestoration by sliding its adjacent fault blocks together. This may joinend points of a horizon H_(τ) separated on opposite sides of the faultin the present day model to merge into the same position in the restoredmodel by sliding the end points towards each other in a directiontangential to the surface of the τ-active fault.

In operation 1660, for each τ-inactive fault, a processor may model theτ-inactive fault, not as a discontinuous fault surface, but as acontinuous non-fault surface in the restoration transformation. Theτ-inactive fault may be modeled as a surface in which the discontinuityinduced by the fault has been deactivated to prevent fault blocks fromsliding in directions tangential to the surface of the fault as timeapproaches the restoration time τ. The processor may model theτ-inactive fault during restoration by keeping collocated points onopposite sides of the fault in the present day model together in therestored model.

After the geological model has been restored for a first pastrestoration time τ (operations 1620-1660), the process may repeat torestore the model for a second different past restoration time τ′. Insome embodiments, the geological model may be sequentially restored to asequence of multiple past restoration times τ₁, τ₂, . . . , τ_(n). Inmultiple (all or not all) of the past restoration times τ₁, τ₂, . . . ,τ_(n), the fault network may be divided into a different subset ofτ-active and τ-inactive faults, e.g., because different faults fracturethe subsurface at different geological times. In some embodiments, aprocessor may play a moving image sequence in which the 3D model isiteratively restored in a forward or reverse order of the sequence ofpast restoration times τ₁, τ₂, . . . , τ_(n) to visualize changes in thesubsurface geology over the passage of time.

In operation 1670, a processor may display a visualization of an imageof the subsurface geology of the Earth overlaid with τ-active faults andτ-inactive faults in the restored model at past restoration time τ. Theprocessor may display the τ-active faults and the τ-inactive faults withdifferent visual identifiers, such as, different levels of translucency,different colors, different patterns, etc.

New Restoration Transformation u_(τ), v_(τ), and t_(τ):

A restoration transformation may transform a geological image of thesubsurface of the Earth from a present day space (e.g., x,y,zcoordinates) to a restoration space (e.g., u_(τ), v_(τ), and t_(τ)coordinates) as it was formed at an intermediate restoration time in thepast τ (before the present-day but after the start of the subsurfacedeposition). An ideal restoration should transform the verticalcoordinate t_(τ) in a manner that strictly honors the thickness oflayers, to preserve areas and volumes of the Earth, so that terrains arenot stretched or squeezed over time in the vertical dimension. However,conventional restoration transformations typically deform the verticalcoordinates, forcing terrains to stretch and squeeze, resulting inerrors in the restoration model.

Embodiments of the invention improve the accuracy of the restorationmodel by establishing a vertical restoration coordinate t_(τ) thatpreserves layer thickness. This may be achieved by implementing athickness-preserving constraint that sets a variation in the verticalrestoration coordinate t_(τ) between any two points separated by aninfinitesimal difference in the direction of maximal variation of thevertical coordinate t_(τ) to be approximately equal to the distancebetween the points in the direction of maximal variation. An example ofthis constraint may be modeled by ∥grad t_(τ)(x,y,z)∥=1. Thisconstraint, however, is non-linear and highly complex and time-consumingto solve. Due to its complexity, this constraint is rarely used inconventional restoration models, and instead replaced byover-simplifications, such as equations (33) and (34), that result inmodel errors as shown in histograms 501 and 502 of FIG. 5, andhistograms 901 and 902 of FIG. 9, respectively.

Embodiments of the invention improve the accuracy of the restored modelby establishing a new thickness-preserving constraint that introduces anerror correction term ε_(τ). The new thickness-preserving constraintsets the restored vertical coordinate t_(τ) to be equal to a sum of afirst approximation t′_(τ) of the vertical coordinate and an errorcorrection term ε_(τ), wherein the error correction term ε_(τ) iscomputed by solving a relationship in which a variation in the sum ofthe first approximation t′_(τ) of the vertical coordinate and the errorcorrection term ε_(τ) between any two points separated by aninfinitesimal difference in the direction of maximal variation of thesum is approximately equal to the distance between the points in thedirection of maximal variation. An example of this constraint may bemodeled by ∥grad (t′_(τ)+ε_(τ))∥=1. The new thickness-preservingconstraint preserves layer thickness with greater accuracy as shown inhistogram 503 of FIG. 5 as compared to conventional approximations shownin histograms 501 and 502 of FIG. 5 and minimizes volume variation withgreater accuracy as shown in histogram 903 of FIG. 9 as compared toconventional approximations shown in histograms 901 and 902 of FIG. 9,respectively.

Embodiments of the invention further improve the performance andcomputational speed of the computer generating the restored model bylinearizing the new thickness-preserving constraint. As an example, thenew thickness-preserving constraint may be linearized as follows. ∥grad(t′_(τ)+ε_(τ))∥=1 may be squared to obtain ∥grad t′_(τ)∥²+∥gradε_(τ)∥²+∥2·grad t′_(τ)·grad ε_(τ)∥=1. The error correction term ε_(τ)may be generated such that the square of its spatial variation, ∥gradε_(τ)∥², is negligible. Accordingly, the thickness-preserving constraintsimplifies to a new linear thickness-preserving constraint of gradε_(τ)·grad t′_(τ)≅½{1−∥grad t′_(τ)∥²} (eqn. (37)). Thisthickness-preserving constraint is linear because t′_(τ) is alreadyknown, so the constraint is a relationship between the gradient of theerror ε_(τ) and the gradient of the known first approximation of thevertical coordinate t′_(τ). The computer may therefore compute the newthickness-preserving constraint in linear time, which is significantlyfaster than computing the non-linear constraints ∥grad t_(τ)∥=1 or ∥grad(t′_(τ)+ε_(τ))∥=1.

Contrary to conventional methods, the computational complexity forperforming the restoration transformation according to embodiments ofthe invention is significantly reduced compared to classical methodsthat are based on the mechanics of continuous media. As a consequence,the modeling computer uses significantly less computational time andstorage space to generate the inventive restoration model.

Contrary to conventional methods that allow variations of geologicalvolumes and deformations, embodiments of the invention implement a newset of geometrical constraints and boundary conditions that preservegeological volumes and deformations while adhering to geologicalboundaries.

Contrary to conventional methods, embodiments of the invention restorefaults along fault striae (e.g., see FIG. 12) induced by the twin pointsassociated with the paleo-geographic coordinates of a depositional(e.g., GeoChron) model, given as input of the restoration method.

An ideal restoration should also transform the horizontal coordinatesu_(τ) and v_(τ) in a manner that strictly honors lateral spatialdistribution, to preserve areas and volumes of the Earth, so thatterrains are not stretched or squeezed over time in the horizontaldimensions. However, conventional restoration transformations based ondepositional coordinates (e.g., paleo-geographic coordinates u and v)typically deform the horizontal coordinates, forcing terrains to stretchand squeeze, resulting in errors in the restoration model.

Embodiments of the invention improve the accuracy of the restorationmodel at time τ by establishing horizontal restoration coordinates u_(τ)and v_(τ) that restore the horizon surface H_(τ) deposited at time τconsistently with horizontal depositional coordinates u and v whilstminimizing deformations. In one embodiment, on the horizon surface H_(τ)only, the horizontal restoration coordinates u_(τ) and v_(τ) are equalto the depositional coordinates u and v (see e.g., equation (20)) andthe spatial variations of the horizontal restoration coordinates u_(τ)and v_(τ) are preserved with respect to the horizontal depositionalcoordinates u and v (see e.g., equation (21)). Thus, each restorationmodel at time τ, presents a horizon surface H_(τ), as it was configuredat that time τ when it was originally deposited. Additionally oralternatively, horizontal restoration coordinates u_(τ) and v_(τ) aremodeled in a tectonic style (e.g., using constraints (22) or (23)) thatis consistent with that of the horizontal coordinates u and v of thedepositional model, which makes the restoration more accurate becausethe geological context is taken into account. Additionally oralternatively, horizontal restoration coordinates u_(τ) and v_(τ) aremodeled to minimize deformations induced by the restoration of horizonH_(τ), rather than minimizing deformations in the whole volume G. Thismay be achieved by implementing constraints (41) and (42) that onlyenforce orthogonality of gradients of u_(τ) and v_(τ) with local axesb_(τ) and a_(τ), but which do not constrain the norm of grad u_(τ) andgrad v_(τ), as is typically constrained for horizontal depositionalcoordinates u and v consistent with the depositional time model.Horizontal restoration coordinates u_(τ) and v_(τ) may also beconstrained only in G_(τ), thereby only taking into account the part ofthe subsurface to be restored, not the entire model G. Additionally oralternatively, horizontal restoration coordinates u_(τ) and v_(τ) may beconstrained to be equal on opposite sides of τ-active faults at twinpoint locations, where the twin points are computed from fault striae,which also ensures consistency with the depositional model (see e.g.,equation (43)). Additionally or alternatively, horizontal restorationcoordinates u_(τ) and v_(τ) are constrained to be equal on oppositesides of τ-inactive faults at mate point locations to cancel the effectof inactive faults on the restoration model (see e.g., equation (43)).

Reference is made to FIG. 17, which is a flowchart of a method torestore a geological model with improved accuracy using a newthickness-preserving constraint, according to an embodiment of theinvention.

In operation 1710, a processor may receive a 3D model of the present daymeasured geometry (e.g., present day model 202) comprising one or morefolded (e.g., curvilinear or non-planar) geological horizon surfaces(e.g., 210). The present day model may be measured tomographically byscanning the Earth's subsurface e.g., as described in reference to FIGS.14 and 15.

In operation 1720, a processor may select or receive a past restorationtime τ that is “intermediate” or prior to the present time and after thestart of the subsurface's deposition (the time period when an oldesthorizon surface in the 3D model was originally deposited).

In operation 1730, a processor may restore the 3D model from the presentday measured geometry (e.g., present day model G_(τ) 202 in xyz-space G220) to the predicted past geometry at the restoration time in the pastτ (e.g., restored model G _(τ) 203 in u_(τ)v_(τ)t_(τ)-space 219) using a3D restoration transformation (e.g., u_(τ)v_(τ)t_(τ)-transform 201). Atthe restored time in the past τ, the geological layers above H_(τ)(e.g., H_(τ+1) . . . H_(n)) did not yet exist, so the subregion aboveH_(τ) in the present day space G 220 is eliminated or omitted, and onlythe subregion G_(τ) 202 below and aligned with H_(τ) (e.g., H₁ . . .H_(τ)) in the present day space G 220 is restored. The 3D restorationtransformation includes a vertical component that restores the geometryto the vertical coordinate t_(τ) and two lateral or horizontalcomponents that restore the geometry to the horizontal coordinates u_(τ)and v_(τ). The restored vertical coordinate t_(τ) and horizontalcoordinates u_(τ) and v_(τ) represent the predicted vertical andhorizontal positions, respectively, where particles in the subsurfacewere located in the Earth at the restoration time in the past τ. Becausethe region of the subsurface deposited after the restoration time τ(e.g., deposited shallower in the subsurface than H_(τ)) did not yetexist at the time of the restored model, the processor may restore andcompute coordinates for the part or subregion G_(τ) of the subsurface Gthat was deposited at a geological time of deposition t prior to orduring the past restoration time τ (e.g., deposited deeper than, or atthe same layer in the subsurface as, H_(τ)). Accordingly, the restoredmodel eliminates or omits all relatively shallower or younger horizonsurfaces or layers that were originally deposited after the pastrestoration time τ.

The processor may restore the vertical coordinate t_(τ) such that pointsalong a horizon surface H_(τ) (e.g., 210) modeling sediment that wasdeposited at the selected restoration time τ have a substantiallyconstant value for the restored vertical coordinate t_(τ) (see e.g.,eqn. (19)). Further, the processor may restore the vertical coordinatet_(τ) such that at any location in the 3D model, the restored verticalcoordinate t_(τ) is equal to a sum of a first approximation t′_(τ) ofthe vertical coordinate and an error correction term ε_(τ), wherein theerror correction term ε_(τ) is computed by solving a relationship inwhich a variation in the sum of the first approximation t′_(τ) of thevertical coordinate and the error correction term ε_(τ) between any twopoints separated by an infinitesimal difference in the direction ofmaximal variation of the sum is approximately equal to the distancebetween the points in the direction of maximal variation. The errorcorrection term ε_(τ) may correct errors in the first approximationt′_(τ) of the vertical coordinate. This constraint may be represented bya linear second order approximation (see e.g., eqn. (37)).

In some embodiments, the processor computes the first approximationt′_(τ) of the vertical coordinate by solving a relationship in which thespatial variation of the vertical coordinate t′_(τ) is locallyapproximately proportional to the spatial variation of a geological timeof deposition t. In some embodiments, the coefficient of proportionalityis locally equal to the inverse of the magnitude of the maximal spatialvariation of the geological time of deposition (see e.g., eqn.(34)-(1)). This relationship may give the vertical restorationcoordinate t_(τ) the shape of the horizon H_(τ) because, on the horizon,the gradient of depositional time t is normal to the horizon surface.Thus, the ratio grad t/∥grad t∥ follows the shape of the horizon.

In some embodiments, the processor computes the first approximationt′_(τ) of the vertical coordinate by solving a relationship in which anyinfinitesimal displacement in the direction orthogonal to horizonsurface H_(τ) results in a variation of the vertical coordinate t′_(τ)approximately equal to the length of the infinitesimal displacement forpoints on the horizon surface H_(τ) (see e.g., eqn. (33)-(1)).

In some embodiments, the processor computes the restored verticalcoordinate t_(τ) in parts of the subsurface which are older thanrestoration time τ such that iso-value surfaces of the restored verticalcoordinate t_(τ) are parallel to the horizon surface H_(τ) and thedifference in the restored vertical coordinate t_(τ) between twoarbitrary iso-values is equal to the distance between the correspondingiso-surfaces (see e.g., eqn. (31)). Parallel surfaces may be planarparallel in the restored model, and curved parallel (e.g., havingparallel tangent surfaces) in present day model, such that the surfacesare non-intersecting at limits.

In some embodiments, the error correction term ε_(τ) is null at pointsalong the horizon surface H_(τ) that was deposited at the selectedrestoration time in the past τ so that the restored horizon surfaceH_(τ) is flat (see e.g., eqn. (36)).

In some embodiments, the restored horizontal coordinates u_(τ) and v_(τ)are constrained such that for each point along the horizon surface H_(τ)that was deposited at the selected restoration time in the past τ: therestored horizontal coordinates u_(τ) and v_(τ) are equal todepositional horizontal coordinates u and v, respectively, and thespatial variations of the restored horizontal coordinates u_(τ) andv_(τ) are equal to the spatial variations of the depositional horizontalcoordinates u and v, respectively (see e.g., eqns. (20)-(21)). Onaverage, globally over the entire model, the processor may compute ∥gradu∥=1 and ∥grad v∥=1. However, locally, this is not necessarily truee.g., on horizon Hτ. So, while the processor sets grad u_(τ)=grad u andgrad v_(τ)=grad v on Hτ, the processor may not constrain ∥grad u_(τ)∥=1and ∥grad v_(τ)∥=1 on Hτ. Moreover, the processor may not constrain gradu_(τ) to be orthogonal to grad t_(τ). This results from the boundarycondition on Hτ and propagation through its constant gradient.

In some embodiments, the restored horizontal coordinates u_(τ) and v_(τ)are constrained in parts of the subsurface which are older thanrestoration time τ such that directions of maximal change of therestored horizontal coordinates u_(τ) and v_(τ) are linearly constrainedby a local co-axis vector b_(τ) and a local axis vector a_(τ),respectively (see e.g., eqn. (41)).

In some embodiments, the local axis vector a_(τ) is orientedapproximately in the direction of maximal change of depositionalhorizontal coordinate u and orthogonal to the direction of maximalchange of the vertical restoration coordinate t_(τ), and the localco-axis vector b_(τ) is oriented orthogonal to the direction of thelocal axis vector a_(τ) and orthogonal to the direction of maximalchange of the vertical restoration coordinate t_(τ) (see e.g., eqn.(40)).

In some embodiments, if the tectonic style of the 3D model is minimaldeformation, the restored horizontal coordinates u_(τ) and v_(τ) arecomputed over the part of the 3D model of the subsurface which is olderthan restoration time τ such that the directions of maximal change ofu_(τ) and v_(τ) are approximately orthogonal to the local co-axis vectorb_(τ) and the local axis vector a_(τ), respectively. For example,equation (40) constrains the local axis vector a_(τ) to be parallel tothe gradient of u and the local co-axis vector b_(τ) to be orthogonal tothe local axis vector a_(τ), which means that the gradient of u isorthogonal to the local co-axis vector b_(τ). Equation (41) furtherconstrains the gradient of u_(τ) to be approximately orthogonal to thelocal co-axis vector b_(τ). Accordingly, the gradient of u_(τ) isapproximately parallel to the gradient of u. The same logic implies thegradient of v_(τ) is approximately parallel to the gradient of v.

In some embodiments, if the tectonic style of the 3D model is flexuralslip, the restored horizontal coordinates u_(τ) and v_(τ) are computedover the part of the 3D model of the subsurface which is older thanrestoration time τ such that projections of their directions of maximalchange over the iso-value surfaces of the restored vertical coordinatet_(τ) are approximately orthogonal to local co-axis vector b_(τ) and thelocal axis vector a_(τ), respectively (see e.g., eqn. (42)).

In some embodiments, the values of the restored horizontal coordinatesu_(τ) and v_(τ) are constrained in parts of the subsurface which areolder than the restoration time τ to be respectively equal on twinpoints on τ-active faults, wherein twin points are points on oppositesides of a τ-active fault that were collocated at the restoration time τand are located on the same fault stria in the present day model, tomerge the twin points into the same position in the restored model bysliding the twin points towards each other in a direction tangential tothe surface of the τ-active fault (see e.g., eqn. (43)).

In some embodiments, the values of the restored horizontal coordinatesu_(τ) and v_(τ) are constrained in parts of the subsurface which areolder than the restoration time τ to be respectively equal on matepoints on τ-inactive faults, wherein mate points are points on oppositesides of a τ-inactive fault that are collocated at present day time, tomove mate points together on opposite sides of τ-inactive faults (seee.g., eqn. (43)).

In operation 1740, a processor may display an image of the restored 3Dmodel of the subsurface geology of the Earth such that each point in the3D model is positioned at the restored coordinates u_(τ), v_(τ), t_(τ)defining the location that a piece of sediment represented by the pointwas located at the restoration time in the past τ.

In some embodiments, the processor may receive an increasingchronological sequence of past restoration times τ₁, τ₂, . . . , τ_(n).For each restoration time τ_(i) in sequence τ₁, τ₂, . . . , τ_(n), theprocessor may repeat operations 1720-1730 to compute a corresponding 3Drestoration transformation Rτ_(i). 3D restoration transformation Rτ_(i)restores the part of the subsurface older than horizon H_(τ) _(i) to itspredicted past geometry at time τ_(i), e.g., to 3D restored coordinatesu_(τ) _(i) , v_(τ) _(i) , and t_(τ) _(i) .

In operation 1750, in some embodiments, a processor may play a movingimage sequence in which the 3D model is iteratively restored in aforward or reverse order of the sequence of past restoration times τ₁,τ₂, . . . , τ_(n) to visualize changes in the subsurface geology overthe passage of time.

In some embodiments, the processor may edit the model in the restorationspace and then reverse the restoration transformation to apply thoseedits in the present day space. For example, the processor may edit thedepositional values u, v, and t associated with the restored 3D model,and then reverse transform the restored 3D model forward in time fromthe predicted past geometry at the restoration time in the past τ to thepresent day measured geometry using an inverse of the 3D restorationtransformation 200 to incorporate the edits from the restored model intothe present day model.

Decompaction at Intermediate Restoration Time τ:

Compaction may refer to the pore space reduction in sediment within theEarth's subsurface. Compaction is typically caused by an increase inload weight of overlying geological layers as they are deposited overtime. As sediment accumulates, compaction typically increases, as timeand depth increase. Conversely, porosity typically decreases, as timeand depth increase. For example, at a depositional time t₀ when a layeris deposited with no overlaying geology, the depositional model hasminimal or no compaction and maximum depositional porosity ψ ₀. At anintermediate restoration time τ, when there is an intermediate load ofoverlying deposited layers, the restored model has an intermediate levelof compaction and an intermediate level of porosity ψ _(τ) (or simplyψ). At the present-day time t_(p), when the present-day model has themost deposited layers, the present-day model typically has a maximallevel of compaction and minimum porosity ψ _(p). Accordingly, thedepositional porosity ψ ₀ is greater than the intermediate time porosityψ _(τ), which in turn is greater than the present-day porosity ψ _(p),i.e., ψ ₀>ψ _(τ)>ψ _(p). Further, because deeper layers are typicallydeposited at relatively earlier times than are shallower layers, withineach model at the same time τ, a relatively deeper geological layertypically experiences a relatively greater load than does a relativelyshallower geological layer, resulting in greater compaction and lesserporosity.

Whereas compaction is a result of deposition over the forward passage oftime, the process of restoration reverses the passage of time tovisualize geology at an intermediate time in the past τ (before thepresent day and after the start of deposition of the oldest subsurfacelayer). Accordingly, embodiments of the invention generate a restorationmodel by reversing the effects of compaction in a process referred to as“decompaction” to more accurately depict how the geometry of geologicallayers change as their depths increase. Whereas compaction compressesthe geological layers, decompaction reverses those effects,decompressing and uplifting terrains, resulting in increased layerthicknesses and increased intermediate time porosity ψ _(τ) (or simplyψ) in the restored domain as compared with the compacted present-daydomain ψ _(p). Decompaction decompresses the geology by a greater amountthe earlier the intermediate restoration time τ is in the past and thedeeper the layer is underneath the Earth's surface.

Conventional decompaction techniques, however, are notoriouslyunreliable. Laboratory experiments on rock samples show that, duringburial when sediments contained in a volume V(r _(τ)) compact undertheir own weight, their porosity ψ(r _(τ)) exponentially decreasesaccording to Athy's law:

ψ(r _(τ))≃ψ ₀(r _(τ))·exp{−κ(r _(τ))·ι(r _(τ))} ∀r _(τ) ∈ G _(τ)  (52)

where V(r _(τ)) represents an infinitely small volume of sedimentcentered on a point r _(τ) ∈ G _(τ) underneath the sea floor Sτ(0) ≡Hτ,δ(r _(τ)) is the absolute distance, or depth, from point r _(τ) ∈ Gτ tosea floor Sτ(0) measured at restoration time τ, and ψ ₀(r _(τ))<1 andκ(r _(r)) are known non-negative coefficients which depend only on rocktype at location r _(τ). ψ ₀(r _(τ)) is the porosity of the rock typewith approximately no (zero) compaction, i.e., the porosity at itsdepositional time t₀ before any layers were deposited to compress fromabove. κ is an experimental measurement derived from compressionexperiments of Athy's law performed in laboratory tests. As an example,assuming that geological depth δ(r _(τ)) is expressed in meters, thefollowing average coefficients for sedimentary terrains were observed insouthern Morocco:

Rock Type Ψ _(o) κ Siltstone 0.62 0.57 × 10⁻³ Clay 0.71 0.77 × 10⁻³Sandstone 0.35 0.60 × 10⁻³ Carbonates 0.46 0.23 × 10⁻³ Dolomites 0.210.61 × 10⁻³Because, in the restored Gτ-space, −t_(τ)(r _(τ)) measures the verticaldistance from point r _(τ) to the sea floor Sτ(0), the depth δ(r _(τ))in equation (52) may be equivalently expressed as:

δ( r _(τ))=−t _(τ)( r _(τ)) ∀ r _(τ) ∈ G _(τ)  (53)

Accordingly, in the context of embodiments of the invention, Athy's lawmay be reformulated as:

$\begin{matrix}{{\overset{\_}{\Psi}\left( {\overset{\_}{r}}_{\tau} \right)} \simeq {{{{\overset{\_}{\Psi}}_{o}\left( {\overset{\_}{r}}_{\tau} \right)} \cdot \exp}\left\{ {{\overset{\_}{\kappa}\left( {\overset{\_}{r}}_{\tau} \right)} \cdot {t_{\tau}\left( {\overset{\_}{r}}_{\tau} \right)}} \right\} {\forall{{\overset{\_}{r}}_{\tau} \in {\overset{\_}{G}}_{\tau}}}}} & (54)\end{matrix}$

Athy's law alone, however, incorrectly models porosity ψ and thereforeoften models decompaction inaccurately. Under Athy's law, restorationporosity ψ depends only on predictions extrapolated based on rockproperties (ψ ₀ and κ), but does not actually measure real-worldporosity. Because Athy's law is not rooted in the real-world geology, itoften leads to inaccurate overestimated or underestimated compaction.Further, Athy’ s law models compaction based on porosity only at thetime of earliest deposition, ψ ₀, but not porosity that occurs in thepresent-day, ψ _(p). Once the model is transformed from the present-dayto restored time τ, but prior to decompaction, the restored model stillerroneously exhibits present-day compaction ψ _(p). Because Athy's lawdoes not eliminate present-day compaction, which erroneouslyover-compresses terrains compared to restoration porosity, the resultingmodel is incorrectly decompacted at the restoration time τ.

Embodiments of the invention improve decompaction techniques by modelingdecompaction at an intermediate restoration time in the past τ based onreal-world measurements of present-day compaction ψ _(p) experimentallyobserved within the subsurface geology of the Earth. Modelingdecompaction based on present-day compaction measurements accounts forthe many real-world geological variables, such as those in the aboveexample scenarios, that Athy's law misses.

Some embodiments accurately decompact the restoration model bysimultaneously (1) removing the impact of present-day compactionaffecting terrains in (incorrectly) restored version at time τ (e.g.,“total” decompaction, such as, defined in equations (58)); and (2)recompacting these terrains according to their depth in the restoredmodel (e.g., “partial” recompaction, such as, defined in equations(59)). Embodiments of the invention solve the difficult problem ofperforming these two operations (decompaction and recompaction)simultaneously.

Reference is made to FIG. 18, which schematically illustrates an example3D geological volume of a compacted model 1810 representing the porosityof a subsurface region before decompaction (right image) and acorresponding decompacted model 1800 representing the porosity of theregion after decompaction (left image) in the restored Gτ-space at anintermediate restoration time in the past τ, according to an embodimentof the invention. Embodiments of the invention replace originalrestoration coordinates of the compacted model 1810 with new restorationcoordinates {u_(τ), v_(τ), t_(τ)}_(r) _(τ) of the decompacted model1800. Decompacted model 1800 may represent a new u_(τ), v_(τ),t_(τ)-transform from the present-day model in Gτ-space to the restoredGτ-space that restores the terrains and induces thickness variations asa consequence of decompaction. This decompaction transformation ismodeled to be the inverse (“reversing time”) of the compaction thatoccurred over the forward passage of time between geological-time τ andthe present geological-time. Some embodiments may start with a region ofthe compacted model 1810 under the horizon Hτ (geology deposited beforetime τ with a present-day level of compaction) and restore the region tothe decompacted model 1800 in Gτ-space (geology deposited before time τwith a level of compaction at intermediate time τ). Because thecompacted model 1810 has not yet been decompacted, its low porosity issimilar to the present-day porosity, yielding vertical lengths dh(r_(τ)) that are too short and compressed for the restored time τ.Accordingly, decompaction vertically stretches the lengths dh⊕(r _(τ))of the decompacted model 1800 to yield a greater porosity predicted atthe time in the past τ. This process may repeat iteratively,layer-by-layer, starting at the top horizon Hτ deposited at therestoration time τ and ending at the bottom horizon deposited at thedepositional time t₀.

Elasto-plastic mechanical frameworks developed to model compaction relyon a number of input parameters which may be difficult for a geologistor geomodeler to assess and are solved using a complex system ofequations. Isostasic approaches are typically simpler to parameterizeand still provide useful information on basin evolution. Therefore,compaction may be considered a primarily one-dimensional verticalcompression induced by gravity which mainly occurs in the early stagesof sediment burial when horizons are still roughly horizontal surfacesclose to the sea floor.

At any point r _(τ) ∈ Gτ within a geological layer, the decompactedthickness dh⊕(r _(τ)) e.g., of a vertical probe of infinitely smallvolume V(r _(τ)) comprising an infinitely short column of sedimentroughly orthogonal to the restored horizon passing through r _(τ) islinked to the thickness dh(r _(τ)) of the shorter, compacted verticalcolumn by, for example, the following relationship:

$\begin{matrix}{{\forall{{\overset{\_}{r}}_{\tau} \in {{\overset{\_}{G}}_{\tau}\text{:}}}}\mspace{14mu} \begin{matrix}{{d{{\overset{\_}{h}}^{\oplus}\left( {\overset{\_}{r}}_{\tau} \right)}} = {{\frac{1}{1 - {{\overset{\_}{\varphi}}_{\tau}\left( {\overset{\_}{r}}_{\tau} \right)}} \cdot d}{\overset{\_}{h}\left( {\overset{\_}{r}}_{\tau} \right)}}} \\{{{with}\text{:}\mspace{14mu} {{\overset{\_}{\varphi}}_{\tau}\left( {\overset{\_}{r}}_{\tau} \right)}} = {{{{\overset{\_}{\Psi}}_{o}\left( {\overset{\_}{r}}_{\tau} \right)} - {{\overset{\_}{\Psi}}_{o}\left( {\overset{\_}{r}}_{\tau} \right)}} \in \left\lbrack {0,1} \right\rbrack}}\end{matrix}} & (55)\end{matrix}$

In this equation, ϕ(r _(τ)) denotes the “compaction coefficient” whichcharacterizes the vertical shortening of the probe at restored locationr _(τ) ∈ Gτ. As an example, FIG. 18 shows the same infinitely shortvertical column of sediment where average porosity is equal to (ψ ₀=⅓)before compaction and (ψ=⅙) after compaction. The compaction coefficient(ψ ₀−ψ) is then equal to (ϕ=⅙) and column shortening (1−ϕ) is (⅚).

Taking Present Day Compaction Into Account to Decompact the RestoredModel in Gτ-Space:

Compacted model 1810, built assuming there is no compaction, incorrectlyignores the compaction characterized by present-day porosity ψ _(p)(r_(τ)). Compacted model 1810 thus results in geology with greatercompaction and smaller porosity than occurred at intermediaterestoration time τ. Embodiments of the invention correct the restoredmodel by decompacting compacted model 1810. The decompaction processinvolves decompressing the vertical dimension's compacted height dh(r_(τ)) or compacted time dt(r _(τ)) (relatively shorter) to elongate thevertical dimension with a decompacted height dh⊕(r _(τ)) or decompactedtime dt⊕(r _(τ)) (relatively longer) (see e.g., equation (60) and/or(64)). This decompaction of height (e.g., in equation (60)) or time(e.g., in equation (64)) is elongated based on compaction coefficient ϕ_(τ) ^(⊖)(r _(τ)), which is a function of the present-day porosity ψ_(p) (r _(τ)) (see e.g., equation (56)). Because the present-dayporosity ψ _(p)(r _(τ)) is less than the restoration porosity ψ _(τ)(r_(τ)), the ratio term in equations (60) and (64) is >1. Accordingly, thedecompacted length dh⊕(r _(τ)) and time dt⊕(r _(τ)) are greater than thecompacted length dh(r _(τ)) and time dt(r _(τ)), respectively, resultingin an elongation of the vertical dimension after decompaction. Thiselongation is thus defined based on real-world measurements of thepresent-day porosity ψ _(p)(r _(τ)), which yields more accuratedecompaction than conventional simulations that ignore real-worldporosity and compaction, such as Athy's law.

Present-day porosity ψ _(p)(r _(τ)) may be measured by direct inspectionof the Earth's subsurface material composition. In one example, porositymay be measured by directly analyzing core samples of the Earth'ssubsurface, for example, using a variety of methods to compare bulk rockvolume and total sample volume. In one example, porosity may be derivedfrom well logs, which are measurements performed on rock inside wells.Samples may be collected and porosity measured at regularly orirregularly spaced intervals within the Earth (e.g., bored into theEarth or along well paths). After porosity measurements are taken atthose discrete locations, porosity may be extrapolated throughout theentire studied domain. In one example, at least one (and preferablymultiple) samples are collected at each distinct depositional layer ordepth (e.g., deposited at each distinct period of time).

Example decompaction processes may proceed as follows:

Let ϕ _(τ) ^(⊖)(r _(τ)) be a total compaction coefficient (representinga total compaction as a difference between the minimum present-dayporosity and maximum depositional porosity) and let ϕ _(τ) ^(⊖)(r _(τ))be an intermediate compaction coefficient (representing a partialcompaction as a difference between the intermediate restoration porosityand maximum depositional porosity). The pair of compaction coefficients,ϕ _(τ) ^(⊖)(r _(τ)) and ϕ _(τ) ^(⊕)(r _(τ)), may be defined, forexample, as:

$\begin{matrix}{{\forall{{\overset{\_}{r}}_{\tau} \in {{\overset{\_}{G}}_{\tau}\text{:}}}}\; \begin{matrix}{{{\overset{\_}{\varphi}}_{\tau}^{\ominus}\left( {\overset{\_}{r}}_{\tau} \right)} = {{{\overset{\_}{\Psi}}_{o}\left( {\overset{\_}{r}}_{\tau} \right)} - {{\overset{\_}{\Psi}}_{p}\left( {\overset{\_}{r}}_{\tau} \right)}}} \\{{{\overset{\_}{\varphi}}_{\tau}^{\oplus}\left( {\overset{\_}{r}}_{\tau} \right)} = {{{\overset{\_}{\Psi}}_{o}\left( {\overset{\_}{r}}_{\tau} \right)} - {\overset{\_}{\Psi}\left( {\overset{\_}{r}}_{\tau} \right)}}}\end{matrix}} & (56)\end{matrix}$

Because compaction typically increases over time, the present dayporosity ψ _(p)(r _(τ)) may be assumed to be less than the restored timeporosity ψ(r _(τ)):

ψ _(p)( r _(τ))<ψ( r _(τ)) ∀ r _(τ) ∈ Gτ  (57)

This inequality implies that intermediate compaction coefficient ϕ _(τ)^(⊕)(r _(τ))≤total compaction coefficient ϕ _(τ) ^(⊖)(r _(τ)), and so,the ratio in equations (60) and (64) is greater than 1, resulting in avertical elongation in height to dh⊕(r _(τ)) and/or time to dt⊕(r _(τ))in the decompacted model 1800 relative to the compacted model 1810.

Considering once again the vertical probe introduced above in restoredspace Gτ, decompaction may proceed by using equation (55) twice, once ina forward and then in a backward transformation, for example, asfollows:

-   -   1. First, to completely cancel out the compaction characterized        by given, present day porosity ψ _(p)(r _(τ)), a “total”        vertical decompaction is applied by updating dh(r _(τ)) as        follows:

$\begin{matrix}{{d{{\overset{\_}{h}}_{o}\left( {\overset{\_}{r}}_{\tau} \right)}} = {{\frac{1}{1 - {{\overset{\_}{\varphi}}_{\tau}^{\ominus}\left( {\overset{\_}{r}}_{\tau} \right)}} \cdot d}{\overset{\_}{h}\left( {\overset{\_}{r}}_{\tau} \right)}}} & (58)\end{matrix}$

-   -   -   After canceling the total compaction in the first operation,            the probe porosity is equal to the depositional porosity ψ            ₀(r _(τ)) having no or negligible compaction.

    -   2. Next, a “partial” recompaction is applied as a function of        the actual restoration porosity ψ(r _(τ)) approximated by        equation (54) at geological-time τ, for example, as:

dh ⊕( r _(τ))={1−ϕ _(τ) ^(⊕)( r _(τ))}·dh ₀( r _(τ))   (59)

-   -   -   After this second operation, the probe porosity is equal to            the intermediate time restoration porosity ψ(r _(τ)).

Therefore, to take present-day compaction into account, equation (55)may be replaced, for example, by:

$\begin{matrix}{{d{{\overset{\_}{h}}^{\oplus}\left( {\overset{\_}{r}}_{\tau} \right)}} = {{\frac{1 - {{\overset{\_}{\varphi}}_{\tau}^{\oplus}\left( {\overset{\_}{r}}_{\tau} \right)}}{1 - {{\overset{\_}{\varphi}}_{\tau}^{\ominus}\left( {\overset{\_}{r}}_{\tau} \right)}} \cdot d}{\overset{\_}{h}\left( {\overset{\_}{r}}_{\tau} \right)}\mspace{14mu} {\forall{{\overset{\_}{r}}_{\tau} \in {\overset{\_}{G}}_{\tau}}}}} & (60)\end{matrix}$

where compaction coefficient ϕ _(τ) ^(⊖)(r _(τ)) is based on themeasured present-day porosity ψ _(p)(r _(τ)), for example, as defined inequation (56). Accordingly, the decompacted vertical thickness dh⊕(r_(τ)) at intermediate restoration time τ is elongated based onreal-world measurements of the present-day compaction ψ _(p)(r _(τ))experimentally observed within the subsurface of the Earth.

Decompaction in GeoChron Based Restoration

In the restored Gτ-space, the geological time of deposition t_(τ)(r_(τ)) may be interpreted as an arc-length abscissa s(r _(τ)) along thevertical straight line passing through r _(τ) oriented in the samedirection as the vertical unit frame vector {r _(t) _(τ) =r _(z)}.Therefore, in the Gτ-space,

dt _(τ)( r _(τ))=ds( r _(τ))=dh ( r _(τ))   (61)

may represent the height of an infinitely short vertical column ofrestored sediment located at point r _(τ)∈Gτ, subject to present-daycompaction. As a consequence, to take compaction into account in therestored Gτ-space, according to equations (60) and (61), geological-timet_(τ)(r _(τ)) may be replaced by a “decompacted” geological-time t_(τ)^(⊕)(r _(τ)) such that, for example:

$\begin{matrix}{{\frac{{dt}_{\tau}^{\ominus}}{{dt}_{\tau}}_{{\overset{\_}{r}}_{\tau}\;}} = {\frac{d{{\overset{\_}{h}}^{\oplus}\left( {\overset{\_}{r}}_{\tau} \right)}}{d{\overset{\_}{h}\left( {\overset{\_}{r}}_{\tau} \right)}} = \frac{1 - {{\overset{\_}{\varphi}}_{\tau}^{\oplus}\left( {\overset{\_}{r}}_{\tau} \right)}}{1 - {{\overset{\_}{\varphi}}_{\tau}^{\ominus}\left( {\overset{\_}{r}}_{\tau} \right)}}}} & (62)\end{matrix}$

Assuming that {r _(t) _(τ) =r _(z)} is the unit vertical frame vector ofthe Gτ-space, it follows, for example, that:

$\begin{matrix}{{{grad}\mspace{14mu} {{t_{\tau}^{\oplus}\left( {\overset{\_}{r}}_{\tau} \right)} \cdot {\overset{\_}{r}}_{t_{\tau}}}} = {{\frac{{dt}_{\tau}^{\ominus}\left( {{\overset{\_}{r}}_{\tau} + {s \cdot {\overset{\_}{r}}_{t_{\tau}}}} \right)}{ds}_{s = 0}} = {\frac{{dt}_{\tau}^{\ominus}}{{dt}_{\tau}}_{{\overset{\_}{r}}_{\tau}\;}}}} & (63)\end{matrix}$

From this, it can be concluded that the compacted geological-timet_(τ)(r _(τ)) of point r _(τ) ∈ Gτ should be transformed into adecompacted geological-time t_(τ) ^(⊕)(r _(τ)), for example, honoringthe following differential equation:

$\begin{matrix}{{{grad}\mspace{14mu} {{t_{\tau}^{\oplus}\left( {\overset{\_}{r}}_{\tau} \right)} \cdot {\overset{\_}{r}}_{t_{\tau}}}} = {\frac{1 - {\varphi_{\tau}^{\oplus}\left( {\overset{\_}{r}}_{\tau} \right)}}{1 - {\varphi_{\tau}^{\ominus}\left( {\overset{\_}{r}}_{\tau} \right)}}\mspace{14mu} {\forall{{\overset{\_}{r}}_{\tau} \in {\overset{\_}{G}}_{\tau}}}}} & (64) \\{{{with}\text{:}\mspace{14mu} {{\overset{\_}{\varphi}}_{\tau}^{\oplus}\left( {\overset{\_}{r}}_{\tau} \right)}} = {{{{{{\overset{\_}{\Psi}}_{o}\left( {\overset{\_}{r}}_{\tau} \right)} - {\overset{\_}{\Psi}\left( {\overset{\_}{r}}_{\tau} \right)}}\&}\mspace{14mu} {{\overset{\_}{\varphi}}_{\tau}^{\oplus}\left( {\overset{\_}{r}}_{\tau} \right)}} = {{{\overset{\_}{\Psi}}_{o}\left( {\overset{\_}{r}}_{\tau} \right)} - {{\overset{\_}{\Psi}}_{p}\left( {\overset{\_}{r}}_{\tau} \right)}}}} & \;\end{matrix}$

Due to the vertical nature of compaction, on the top restored horizon{S(0)≡H _(τ)}, geological-time t_(τ) ^(⊕)(r _(τ)) should vanish orreduce to zero and its gradient should be vertical. In other words, inaddition to the constraint of equation (64), geological-time t_(τ)^(⊕)(r _(τ)) may also honor the following example boundary conditionswhere r _(u) _(r) and r _(v) _(r) may represent the unit horizontalframe vectors of the Gτ-space:

$\begin{matrix}{{\forall{{\overset{\_}{r}}_{\tau}^{o} \in {\left\{ {{{\overset{\_}{}}_{\tau}(0)} \equiv {\overset{\_}{H}}_{\tau}} \right\} \text{:}}}}\mspace{14mu} \begin{matrix}\left. 1 \right) & {{t_{\tau}^{\oplus}\left( {\overset{\_}{r}}_{\tau}^{o} \right)} = 0} \\\left. 2 \right) & {{{grad}\mspace{14mu} {{t_{\tau}^{\oplus}\left( {\overset{\_}{r}}_{\tau}^{o} \right)} \cdot {\overset{\_}{r}}_{{\overset{\_}{u}}_{\tau}}}} = 0} \\\left. 3 \right) & {{{grad}\mspace{14mu} {{t_{\tau}^{\oplus}\left( {\overset{\_}{r}}_{\tau}^{o} \right)} \cdot {\overset{\_}{r}}_{{\overset{\_}{v}}_{\tau}}}} = 0}\end{matrix}} & (65)\end{matrix}$

Boundary condition (65)(1) may ensure that the top restored horizon H_(τ) is flat and planar at intermediate restoration time τ when it wasdeposited. Boundary conditions (65)(2) and (65)(3) may ensure that thedirection of change (gradient) of the geological-time t_(τ)(r _(τ)) isvertical in the Gτ-space.

As compaction is a continuous process, geological-time t_(τ) ^(⊕)(r_(τ)) may be continuous (e.g., C⁰-continuous) across all faultsaffecting Gτ. As a consequence, in addition to the constraints inequations (64) and (65), for any fault F in Gτ, geological-time t_(τ)^(⊕)(r _(τ)) may also honor the following boundary conditions where (r_(F) ^(⊕), r _(F) ^(⊖)) are pairs of “τ-mate-points” defined ascollocated points respectively lying on the positive face F ⁺ andnegative face F ⁻ (opposite sides of fault F) at geological time τ:

t _(τ)( r _(F) ^(⊕))=t _(τ)( r _(F) ^(⊖))

∀F ∈ G _(τ) & ∀(r _(F) ^(⊕), r _(F) ^(⊖))_(τ) ∈ F  (66)

Boundary condition (66) may ensure that, for any pair of collocatedpoints on opposite sides of the fault, the two points have the samedecompacted geological-time coordinate t_(τ) ^(⊕)(r _(τ)). This ensuresthere are no (or reduced) gaps or overlaps along the fault in therestored Gτ-space.

Using an appropriate numerical method, t_(τ) ^(⊕)(r _(τ)) may becomputed in Gτ whilst ensuring that differential equation (64) andboundary conditions (65) and (66) are honored. To ensure smoothness anduniqueness of t_(τ) ^(⊕)(r _(τ)), the following constraint may also beadded:

$\begin{matrix}{\sum\limits_{{({a,b})} \in {\{{u_{\tau},v_{\tau},t_{\tau}}\}}^{3}}{\int_{{\overset{\_}{G}}_{\tau}}^{\;}{{\left\{ {\partial_{a}{\partial_{b}{t_{\tau}^{\oplus}\left( {\overset{\_}{r}}_{\tau} \right)}}} \right\}^{2} \cdot d}{\overset{\_}{r}}_{\tau}\mspace{14mu} {minimum}}}} & (67)\end{matrix}$

In summary, the following GeoChron Based Restoration technique may beused to take compaction into account:

-   -   1. Compute a numerical approximation of the elongated        geological-time t_(τ) ^(⊕)(r _(τ)) in Gτ and use the reverse        u_(τ), v_(τ), t_(τ)-transform to update t_(τ)(r_(τ)) in Gτ:

t _(τ)(r _(τ))←t _(τ) ^(⊖)( r _(τ)) ∀ r _(τ) ∈ G _(τ);   (68)

-   -   2. Recompute numerical approximations of restoration functions        u_(τ)(r_(τ)) and v_(τ)(r_(τ)) in Gτ to prevent voids and        overlaps from being generated in the restored space, as,        according to equations (22) and (23), u_(τ)(r_(τ)) and        v_(τ)(r_(τ)) depend on t_(τ)(r_(τ));    -   3. Build the “decompacted” restored space Gτ as the new, direct        u_(τ), v_(τ), t_(τ)-transform of geological space Gτ observed        today.

This approach to decompaction may be seamlessly integrated into theGeoChron Based Restoration framework according to embodiments of theinvention and is wholly dissimilar to the sequential decompactionfollowing Athy's law along IPG-lines. In particular, embodiments of theinvention perform decompaction based on real-world present-day porosity,a quantity that is accurately measured and extrapolated for any type ofrock without having to make assumptions. Additionally, embodiments ofthe invention allow decompaction in the restored Gτ space representingthe Earth's subsurface at an intermediate restoration time in the pastτ, before the present day and after the start of deposition of theoldest subsurface layer being imaged.

An Analytical Solution

In the general case, the system of equations (64), (65) and (66) istypically too complex to be solved analytically and may be approximatedusing numerical methods. However, in a specific case where κ(r _(τ)), ψ₀(r _(τ)) and ψ _(p)(r _(τ)) are all constant ∀ r _(τ) ∈ Gτ, thecompaction ratio may be integrated at once over the entire domain, andthere is no need to iteratively and independently decompact one layer ata time. This special case allows an analytical solution to the system ofequations, for example, as follows.

In this special case, in Gτ, terrain porosity is homogeneous andcharacterized as for example follows where κ, ψ₀ and ψ_(p) are knownconstants:

$\begin{matrix}{{\forall{{\overset{\_}{r}}_{\tau} \in {{\overset{\_}{G}}_{\tau}\text{:}}}}\mspace{14mu} \begin{matrix}{{\overset{\_}{\kappa}\left( {\overset{\_}{r}}_{\tau} \right)} = \kappa} \\{{{\overset{\_}{\Psi}}_{o}\left( {\overset{\_}{r}}_{\tau} \right)} = \Psi_{o}} \\{{{\overset{\_}{\Psi}}_{p}\left( {\overset{\_}{r}}_{\tau} \right)} = \Psi_{p}}\end{matrix}} & (69)\end{matrix}$

Due to its homogeneity, Gτ may be considered continuous and theintrinsic, vertical nature of compaction implies that any function φ(r_(τ)) defined in Gτ associated to compaction may only depend on thevertical component t_(τ)(r _(τ)) of r _(τ). Therefore, it follows, forexample, that:

φ( r _(τ))=φ(t _(τ)( r _(τ)))=φ(t _(τ))   (70)

Let constants A and B be defined, for example, by:

$\begin{matrix}{{A = \frac{1 - \Psi_{o}}{1 - \Psi_{o} + \Psi_{p}}};{B = \frac{\Psi_{o}}{1 - \Psi_{o} + \Psi_{p}}}} & (71)\end{matrix}$

Let the following example functions be derived from Athy's law inequation (54) and equations (56) and (69):

ϕ _(τ) ^(⊖)(t _(τ))=Ψ_(o)−Ψ_(p); ϕ _(τ) ^(⊕)(t _(τ))=Ψ_(o)−Ψ_(p)·exp(κ·t_(τ)) ∀ t _(τ)≤0   (72)

On the one hand, the following example indefinite integral holds true:

$\begin{matrix}{{\int{\frac{1 - {{\overset{\_}{\varphi}}_{r}^{\oplus}(s)}}{1 - {{\overset{\_}{\varphi}}_{r}^{\ominus}(s)}} \cdot {ds}}} = {{\int{\left\{ {A + {B \cdot {\exp \left( {\kappa \; s} \right)}}} \right\} \cdot {ds}}} = {{A \cdot s} + {\frac{B}{\kappa} \cdot {\exp \left( {\kappa \; s} \right)}}}}} & (73)\end{matrix}$

On the other hand, according to equation (62):

$\begin{matrix}{\; {{t_{\tau}^{\oplus}\left( t_{\tau} \right)} = {{\int_{0}^{t_{\tau}}{\frac{1 - {{\overset{\_}{\varphi}}_{r}^{\oplus}(s)}}{1 - {{\overset{\_}{\varphi}}_{r}^{\ominus}(s)}} \cdot {ds}}} = {{A \cdot t_{\tau}} + {\frac{B}{\kappa} \cdot \left\{ {{\exp \left( {k \cdot t_{\tau}} \right)} - 1} \right\}}}}}} & (74)\end{matrix}$

Therefore, for any {t_(τ)≤0}, the decompacted restoration function t_(τ)^(⊕)(t_(τ)) may be analytically defined, for example, by:

$\begin{matrix}{\; {{t_{\tau}^{\oplus}\left( t_{\tau} \right)} = {{A \cdot t_{\tau}} + {\frac{B}{\kappa} \cdot \left\{ {{\exp \left( {k \cdot t_{\tau}} \right)} - 1} \right\}}}}} & (75) \\{{{{with}\text{:}\mspace{14mu} A} = \frac{1 - \Psi_{o}}{1 - \Psi_{o} + \Psi_{p}}};{B = \frac{\Psi_{o}}{1 - \Psi_{o} + \Psi_{p}}}} & \;\end{matrix}$

Other equations or permutations of these equations or terms may also beused.

Reference is made to FIG. 19, which is a flowchart of a method fordecompacting a 3D model of the subsurface geology of the Earth at anintermediate restoration time in the past τ, according to an embodimentof the invention.

In operation 1910, a processor may receive a 3D model of present-daygeometry of the subsurface geology and a measure of present-day porosityexperimentally measured within the subsurface geology of the Earth. Thepresent day model may be measured tomographically by scanning theEarth's subsurface e.g., as described in reference to FIGS. 14 and 15.To obtain the measure of present-day porosity, a probe may burrow intothe Earth's subsurface or into one or more wells to collect and/oranalyze material from within the subsurface geology of the Earth.Samples of subsurface materials are collected at spaced intervals, fromwhich porosity is extrapolated throughout the studied domain.

In operation 1920, a processor may select or receive a past restorationtime τ that is “intermediate” or prior to the present time and after thestart of the subsurface's deposition (the time period when an oldesthorizon surface in the 3D model was originally deposited).

In operation 1930, a processor may restore the 3D model from the presentday measured geometry (e.g., present day model G_(τ) 202 in xyz-space G220) to the predicted past geometry at the restoration time in the pastτ (e.g., restored model G _(τ) 203 in u_(τ)v_(τ)t_(τ)-space 219) using a3D restoration transformation (e.g., u_(τ)v_(τ)t_(τ)-transform 201). The3D model may be restored, for example, as described in reference to FIG.17. Prior to decompaction, the restored model may be a compacted model(e.g., 1810 of FIG. 18).

In operation 1940, a processor may decompact the vertical dimension ofthe restored 3D model. This may expand, stretch or elongate compactedvertical lengths in the compacted model (e.g., 1810 of FIG. 18) torelatively longer vertical lengths in a decompacted model (e.g., 1800 ofFIG. 18). In one embodiment, the vertical length may be a measure ofheight and the vertical dimension may be expanded from relativelyshorter heights dh(r _(τ)) in the compacted model to relatively longerheights dh⊕(r _(τ)) in a decompacted model (e.g., as defined in equation(60)). In another embodiment, the vertical length may be a measure ofgeological-time when the particles of sediment were originally depositedon the Earth's surface and the vertical dimension may be expanded fromrelatively shorter times dt(r _(τ)) in the compacted model to relativelylonger times dt⊕(r _(τ)) in a decompacted model (e.g., as defined inequation (64)). The vertical lengths may be elongated based on arelationship between a depositional porosity (e.g., ψ ₀(r _(τ))) of thegeological layers at the time sediment in those layers was deposited,restoration porosity (e.g., ψ(r _(τ))) of the geological layers at therestoration time in the past τ, and the present-day porosity (e.g., ψ_(p)(r _(τ))) of the geological layers experimentally measured in thepresent-day. In some embodiments, the relationship between thedepositional porosity, the restoration porosity, and the present-dayporosity may be, for example:

$\frac{1 - {\overset{\_}{\varphi}}_{\tau}^{\oplus}}{1 - {\overset{\_}{\varphi}}_{\tau}^{\ominus}},$

where compaction coefficients ϕ _(τ) ^(⊕)=ψ ₀−ψ and ϕ _(τ) ^(⊖)=ψ ₀−ψ_(p) for all points in the restored 3D model, e.g., as defined inequations (56), (60) and (64). Since porosity decreases over time, therestoration porosity is typically greater than the present-day porosity(e.g., equation (57)) and typically less than the depositional porosity.Accordingly, the compaction coefficients have a relationship ϕ _(τ)^(⊕)≤ϕ _(τ) ^(⊖), and the relationship between the depositional,restoration, and present-day porosities, e.g.,

$\frac{1 - {\overset{¯}{\varphi}}_{\tau}^{\oplus}}{1 - {\overset{\_}{\varphi}}_{\tau}^{\ominus}},$

is greater than 1, resulting in a stretching or elongating effect toincrease the vertical lengths when they are decompacted.

In some embodiments, a processor may decompact the vertical dimension ofthe restored 3D model by a combination (e.g., equation (60)) of totaldecompaction corresponding to an increase in porosity from the presentday porosity to the depositional porosity (e.g., equation (58)) andpartial recompaction corresponding to a partial decrease in the porosityfrom the depositional porosity to the restored porosity (e.g., equation(59)).

At the restored time in the past z, the geological layers above H _(τ)(e.g., H _(τ+1) . . . H _(n)) did not yet exist, so decompaction mayelongate lengths of geological layers below the horizon layer H _(τ)(e.g., H ₀ . . . H _(τ)) deposited at the restoration time in the pastτ. In some embodiments, decompaction may be performed by iterativelydecompacting the subsurface layer-by-layer, starting at the top horizonH _(τ) deposited at the restoration time τ and ending at the bottomhorizon H ₀ deposited at the depositional time. In some embodiments, thedepositional porosity and the present-day porosity may be independentlydetermined for each geological layer of the subsurface. In otherembodiments, when the depositional porosity and the present-day porosityare substantially constant throughout the subsurface geology,decompaction may occur in one operation over the entire domain of therestored 3D model (e.g., as in equation (76)).

Some embodiments may implement a boundary condition that ensures that atop horizon H _(τ) deposited at the restoration time τ is a horizontalplane in the restored 3D model (e.g., equation (65)(1)). Additionally oralternatively, some embodiments may implement a boundary condition thatensures that a direction of change of geological-time when the particlesof sediment were originally deposited on the Earth's surface is verticalin the restored 3D model (e.g., equation (65)(2) and (65)(3)).Additionally or alternatively, some embodiments may implement a boundarycondition that ensures that, for any pair of collocated points onopposite sides of a fault, the two collocated points are decompacted tohave the same coordinate (e.g., equation (66)).

In some embodiments, for example, implemented in a past-time model, suchas the GeoChron model, a processor may decompact the vertical dimensionof the restored 3D model by: computing an elongated geological-time(e.g., dt⊕(r _(τ))) in the restored 3D model (e.g., by solving equation(64)), transforming the elongated geological-time from the restored 3Dmodel to generate a 1D geological-time (e.g., t_(τ)(r_(τ))) in thepresent-day 3D model (e.g., equation (68)), computing 2Dpaleo-depositional coordinates (e.g., u_(τ)(r_(τ)) and v_(τ)(r_(τ)))based on the transformed geological-time (e.g., t_(τ)(r_(τ))) in thepresent-day 3D model, and performing a 3D transformation (e.g., au_(τ),v_(τ),t_(τ)-transformation)) comprising the 1D geological-time and2D paleo-depositional coordinates from the present-day 3D model (e.g.,G_(τ)) to the restored 3D model (e.g., G _(τ)) that is decompacted basedon the elongated geological-time.

Operations of FIGS. 16, 17 and 19 may be performed for example usingsystem 1505 of FIG. 15, e.g., by one or more processor(s) 140 of FIG.15, or another suitable computing system. The embodiments disclosed inreference to FIGS. 16, 17 and 19 may be performed using other operationsor orders of the operations, and the exact set of steps shown in thefigures may be varied.

In the past 30 years, many methods have been proposed to buildgeological models of sedimentary terrains having layers that are bothfolded and faulted. For any given geological-time τ, checking geologicalmodel consistency is considered both simpler and more accurate ifterrains have previously been “restored” to their pre-deformational,unfolded and unfaulted state, as they were at geological-time τ.

Embodiments of the invention provide a new, purely geometrical 3Drestoration method based on the input of a depositional (e.g., GeoChronmodel). Embodiments of the invention are able to handle depositionalmodels of any degree of geometrical and topological complexity, withboth small and large deformations, do not assume elastic mechanicalbehavior, and do not require any prior knowledge of geo-mechanicalproperties. Embodiments of the invention further reduce or eliminategaps and overlaps along faults as part of the restoration transformationand do not resort to any post-processing to minimize such gaps andoverlaps. Compared to other conventional methods, embodiments of theinvention minimize deformations and volume variations induced bygeological restoration with a higher degree of precision, unequaled sofar (see e.g., FIG. 5 and FIG. 9). Embodiments of the invention furtherensure that 2D deformations of horizon surfaces induced by theuvt-transform are kept coherent with 3D deformations of volumes inducedby the new proposed 3D restoration method.

Referring to FIG. 3, for a given restoration time τ, the set of faultsis split into τ-active and τ-inactive subsets. Such a distinctionallows:

-   -   deformations along faults 105 to be minimized,    -   restoration to work correctly even though there are regions of        G_(τ) not continuously connected to H_(τ),    -   gaps and overlaps along faults and the geometry of fault striae        600 are minimized by the restoration transformation, so no        post-processing is needed to correct gaps or overlaps.

Embodiments of the invention input a 3D model of sedimentary terrains inthe subsurface. In one example, the input model may be the GeoChron™model generated by SKUA® software for use in mining and oil and gasindustries. Embodiments of the invention may build a 3D restorationtransformation of this model in such a way that, after transformation,the new model represents terrains as they were at a given intermediaterestoration-time τ (where τ₁<τ<τ₂, before the present day τ₂ and afterthe time of the deposition of the oldest layer τ₁).

For example, G may represent the present day 3D geological domain of theregion of the subsurface being modeled and G_(τ) 202 may represent thesubset of G containing particles of sediment that were deposited at atime prior to or equal to τ. In some embodiments, for all points r ∈ G,a geologic restoration transformation may move a particle of sedimentobserved today at location r to a new restored location r _(τ)(r), e.g.,defined as follows:

r _(τ)(r)=r+R _(τ)(r) ∀ r ∈ G _(τ)  (1)

where R_(τ)(r) represents a 3D field of restoration vectors, e.g.,generated to minimize deformations in G_(τ).

Depositional Model

A depositional model may be generated by inputting a tomographic modelof the present day subsurface geology of the Earth and transforming thatgeology to a past depositional time as each particle was configured whenoriginally deposited in the Earth. Sedimentary particles are depositedover time in layers from deepest to shallowest from the earliest to themost recent geological time periods. Since various layers of terrain aredeposited at different geological times, a depositional model does notimage the geology at any one particular time period, but across manytimes periods, each layer modeled at the geological time when the layerwas deposited. Accordingly, the vertical axis or depth in thedepositional model may be a time dimension representing the time periodof deposition, progressing from oldest to newest geological time as themodel progresses vertically from deepest to shallowest layers.

In one embodiment, the depositional model may be the GeoChron™ model,which is generated by SKUA™ software, that is routinely used by many oil& gas companies to build models of geologic reservoirs which helpoptimize hydrocarbon production and exploration. An exampleimplementation of the GeoChron model is discussed in U.S. Pat. No.8,600,708, which is incorporated by reference herein in its entirety.The depositional model is described in reference to the GeoChron modelonly for example, though any other depositional model may be used.

Reference is made to FIG. 7, which schematically illustrates an exampletransformation from a present day model (upper-left image) to adepositional GeoChron model (bottom-right image), according to anembodiment of the invention. The transformation may be referred to as a“uvt-transform” 700 that transforms a particle of sediment observedtoday at location r=r(x,y,z) in the present day geological domain G(also referred to as “G-space”) 220 to be moved to a new depositionallocation r(r)=r(u,v,t) in the depositional geological domain G (alsoreferred to as “G-space”). The new depositional location r has avertical coordinate that is the geological time t(r) when the particleat location r was deposited and has horizontal or paleo-geographiccoordinates {u(r), v(r)} equal to the lateral spatial location where theparticle at r was located at its depositional time t(r). Thepaleo-geographic coordinates {u(r), v(r)} may be linked to the verticaltime coordinate t(r) by different relationships (e.g., constrained bydifferent systems of differential equations) depending on the structuralstyle of their deposition (e.g., minimal deformation or flexural slip).

In the example uvt-transform 700 shown in FIG. 7, when the geologicaltime coordinate t(r) is equal to the curvilinear distance to the tophorizon H_(τ) 210 along curvilinear axis 222, the uvt-transform is avalid technique for imaging the depositional model. In other words, theuvt-transform is a valid depositional rendering technique if the moduleof its gradient grad t(r) honors the following constraint:

∥grad t(r)∥=1 ∀ r ∈ G   (3)

Embodiments of the invention observe that when ∥grad t(r)∥ differs from“1,” replacing the depositional coordinates {u(r), v(r), t(r)} of theuvt-transform 700 by new restoration coordinates {u_(τ)(r), v_(τ)(r),t_(τ)(r)} where ∥grad t_(τ)∥=1 allows the uvt-transform to be replacedby a u_(τ)v_(τ)t_(τ)-transform that generates a valid restoration modelat restoration time_(τ).

In some embodiments, the depositional (e.g., GeoChron) model includesthe following data structures stored in a memory (e.g., memory 150 ofFIG. 15) (see FIGS. 1, 3, 6, and 7):

-   -   A network of geological faults 105 within the present day domain        G 220.    -   A 3D corner-point grid Γ 100 that fills the G-space 220 with 3D        polyhedral cells 108 (e.g., tetrahedra or hexahedra), without        any gaps or overlaps in the studied domain, in such a way that        no cell edge crosses any fault. The location of each node α 107        of grid Γ 100 in the G-space is denoted r(α).    -   For each geological fault F 105, two disconnected and        independently meshed, collocated surfaces F⁺ 103 and F⁻ 104 on        opposite sides of the fault 105. Surfaces F⁺ 103 and F⁻ 104 may        be composed of 2D facets from the 3D polyhedral cells of grid Γ        100 bordering F 105 on either side of the fault 105. Fault        surfaces F⁺ 103 and F⁻ 104 that are collocated in the present        day model may, during the restoration process of transforming        the model backwards in time, typically slide against one        another, without generating gaps or overlaps between adjacent        fault surfaces or fault blocks.    -   Referring to FIGS. 1, 3 and 6, for each fault F 105, a set of        pairs of points (r_(F) ⁺, r_(F) ⁻) (101, 102) called        “twin-points,” such that:        -   1. The two twin points in each pair are located on opposite            sides of a corresponding pair of twin fault surfaces, r_(F)            ⁺ ∈ F⁺ and r_(F) ⁻ ∈F⁻.        -   2. At geological times before fault F 105 formed in the            subsurface, particles of sediment were collocated which are            observed today at locations r_(F) ⁺ and r_(F) ⁻.    -   During the activation of fault F 105, particles of sediment        initially located on F are assumed to slide along fault-striae        (e.g., see FIG. 12), which are the shortest paths, on F, between        pairs of twin points (r_(F) ⁺, r_(F) ⁻) (101, 102).    -   A tectonic style which may be either a “minimal deformation”        style or a “flexural slip” style. Choosing this tectonic style        is a model decision assumed to have been made by a structural        geologist.    -   A triplet (e.g., {u(r), v(r), t(r)}) of discrete coordinates        defined on a 3D grid Γ 100 of the depositional G-space, such        that, for a particle of sediment observed today at location r,        the coordinate values {u(r), v(r)} represent the        paleo-geographic coordinates of the particle at geological-time        t(r) during the time period when it was deposited. According to        the depositional (e.g., GeoChron) model, the paleo-geographic        coordinates {u(r), v(r)} may honor different differential        equations depending on the tectonic style.

Moreover, referring to FIG. 1 and FIG. 2, the depositional model mayhave the following properties:

-   -   Within the present day domain G, each geological horizon H_(τ)        210 may be defined by a set of particles of sediment which were        deposited at geological time τ:

r ∈ H _(τ) ⇐⇒ t(r)=τ  (4)

-   -   In other words, each horizon H_(τ) 210 is a level-set (constant        value) surface of the geological-time t.    -   Paleo-geographic coordinates {u(r), v(r)} and twin-points (101,        102) given as input are linked e.g. by the following equations:

$\begin{matrix}\left\{ {\left( {r_{F}^{+},r_{F}^{-}} \right)\mspace{14mu} \text{is a pair of twin-points\}}}\; \right. \\\left. \Leftrightarrow\left\{ {\begin{matrix}\begin{matrix}\begin{matrix}{{{{r_{F}^{+} \in F^{+}}\mspace{14mu}\&}\mspace{31mu} r_{F}^{-}} \in F^{-}} \\{{u\left( r_{F}^{-} \right)} = {u\left( r_{F}^{+} \right)}}\end{matrix} \\{{v\left( r_{F}^{-} \right)} = {v\left( r_{F}^{+} \right)}}\end{matrix} \\{{t\left( r_{F}^{-} \right)} = {t\left( r_{F}^{+} \right)}}\end{matrix}\mspace{371mu} \begin{matrix}(5) \\(6) \\(7) \\(8)\end{matrix}} \right. \right.\end{matrix}$

-   -   Additionally or alternatively, each pair of twin-points (r_(F)        ⁺, r_(F) ⁻) (101, 102) may be the intersection of a level set        210 of vertical depositional coordinate t(r) with a “fault        stria” σ(r_(F) ⁻) 600 comprising a curved surface passing        through point r_(F) ⁻ 102 whose geometry is defined by        geological rules, e.g., defining fault blocks sliding against        one another according to tectonic forces and geological context.        As a consequence of constraints defined by equations (6), (7),        and (8) above, fault-striae (e.g., see FIG. 12) may characterize        the paleo-geographic coordinates {u(r), v(r)}, and vice versa.    -   Each point r ∈ G 214 may be characterized by its present day        coordinates (e.g., {x(r), y(r), z(r)}) with respect to a present        day coordinate system {r_(x), r_(y), r_(z)} 220 comprising three        mutually orthogonal unit vectors, e.g., where r_(z) is assumed        to be oriented upward.

It would be appreciated by a person of ordinary skill in the art thatthe GeoChron model and its features described herein are discussed onlyas an example of a depositional model and that these elements may differin other models or implementations without changing the essence of theinvention.

u_(τ)v_(τ)t_(τ)-Transformation

Referring to the volume deformation of FIG. 2, the restoration time τmay be a given geological time in the past and subdomain G_(τ) 202 maybe a part of a 3D present day geological domain G that has terrainsolder than (deposited at a time prior to) or equal to restoration time τand defined by a depositional model. Embodiments of the inventionprovide a new unfolding technique that replaces the uvt-transform(converting the present day model to a depositional model rendering alllayers at their many respective times of deposition) by au_(τ)v_(τ)t_(τ)-transform 201 (converting the present day model to arestored model at a single restoration time τ before the present day butafter the earliest times of deposition of the deepest model layer):

$\begin{matrix}{\left( {x,y,z} \right)\overset{u_{\tau}v_{\tau}t_{\tau}}{\rightarrow}\left\{ {{u_{\tau}\left( {x,y,z} \right)},\ {v_{\tau}\left( {x,y,z} \right)},\ {t_{\tau}\left( {x,y,z} \right)}} \right\}} & (9)\end{matrix}$

Accordingly, present day geological space G_(τ) 202 is transformed intoa restored geological space G _(τ) 203, such that:

-   -   t_(τ)(r) is a vertical spatial coordinate of the subsurface at        the past restoration time τ, and is derived from, but different        than, the geological time of deposition t(r). The vertical        restoration coordinate t_(τ)(r) honors the following constraint:

∥grad t _(τ)(r)∥=1 ∀r ∈ G _(τ)  (10)

-   -   {u_(τ)(r), v_(τ)(r)} are lateral restoration coordinates derived        from, but different than, the paleogeographic coordinates {u(r),        v(r)} of the depositional model.    -   restoration coordinates {u_(τ)(r), v_(τ)(r), t_(τ)(r)} honor        specific inventive constraints described below in such a way        that, using the u_(τ)v_(τ)t_(τ)-transform as a restoration        operator minimizes deformations in the present day domain G_(τ).    -   for each point r ∈ G_(τ) 202, the restoration vector field R_(τ)        may be defined e.g. by:

R _(τ)(r)=u _(τ)(r)·r _(x) +v _(τ)(r)·r _(y) +t _(τ)(r)·r _(z) −r   (11)

Volume Deformation

Compaction may be handled in pre and post-restoration stages, as isknown in the art. Thus, the model may be restored without takingcompaction into account.

Some embodiments of the invention provide an inventive volumedeformation with a new set of inventive geometric constraints on thedepositional model to allow geologic layers to be restored at a givengeological time τ with a precision that has never before been reached.As shown in FIG. 2, in this volume deformation, paleo-geographiccoordinates {u(r) and v(r)} and the geological time coordinate t(r) arereplaced by new restoration coordinates respectively denoted {u_(τ)(r),v_(τ)(r)} and t_(τ)(r).

As shown in FIG. 2:

-   -   a restored volume 203, denoted G _(τ) contains a (e.g., direct        or “right-handed”) coordinate space 219 having orthogonal        coordinate unit vectors {r _(u) _(τ) , r _(v) _(τ) , r _(t) _(τ)        } and a family of horizontal planes {S _(τ)(d):d≥0} 207 parallel        to horizontal coordinate vectors {r _(u) _(τ) , r _(v) _(τ) };    -   a deformed version G_(τ) of G _(τ) contains a (e.g., direct or        “right-handed”) coordinate space 220 having orthogonal        coordinate unit vectors {r_(x), τ_(y), τ_(z)} and a family of        curved surfaces {S_(τ)(d):d≥0} 208 parallel to horizon        {H_(τ)≡S_(τ)(0)} (210=208).

For simplicity and without loss of generality, the coordinate frame unitvectors {r _(u) _(τ) , r _(v) _(τ) , r _(t) _(τ) } 219 of the G_(τ)-space and its origin Ō_(u) _(τ) _(v) _(τ) _(t) _(τ) may be equal tothe coordinate frame unit vectors {r_(x), r_(y), r_(z)} 220 of theG-space and its origin O_(xyz):

r _(u) _(τ) ≡r_(x); r _(v) _(τ) ≡r_(y); r _(t) _(τ) ≡r_(z)

Ō_(u) _(τ) _(v) _(τ) _(t) _(τ) ≡O_(xyz)   (12)

Referring to FIG. 2, the following notation is used:

-   -   each point r ∈ G_(τ) 214 is transformed into point r _(τ) ∈ G        _(τ) 213 and vice versa:

r _(τ) ∈ G _(τ)←→r ∈ G_(τ)  (13)

-   -   {u_(τ)(r _(τ)), v_(τ)(r _(τ))} represent the horizontal        restoration coordinates of r _(τ) 213 with respect to {r _(u)        _(τ) , r _(v) _(τ) } 219 whilst t_(τ)(r) represents the vertical        restoration coordinate or altitude 204 of r _(τ) 213 with        respect to r _(t), oriented upward;    -   for each point r ∈ G_(τ) 214:        -   {x(r), y(r)} represent the horizontal present day            coordinates of r with respect to {r_(x), τ_(y)} whilst z(r)            represents the vertical present day coordinate or altitude            of r with respect to the vertical unit frame vector r_(z)            oriented upward; and        -   {u_(τ)(r), v_(τ)(r), t_(τ)(r)} represent the restoration            coordinates of r _(τ) with respect to the restoration            coordinate frame unit vectors {r _(u) _(τ) , r _(v) _(τ) , r            _(t) _(τ) } 219 of the restored volume G _(τ) 203.

Equivalently to equations (12) and in accordance with equation (1),during restoration of G_(τ), a particle of sediment observed today atlocation r 214 is moved to a new location {circumflex over (r)}(r) 213defined e.g., as follows, where R_(τ)(r) is a restoration vector field:

r (r)=r+R _(τ)(r)   (14)

with, in matrix notation:

$\begin{matrix}{{R_{\tau}(r)} = {\left\lbrack {r_{x},r_{y},r_{z}} \right\rbrack \cdot \begin{bmatrix}{{u_{\tau}(r)} - {x(r)}} \\{{v_{\tau}(r)} - {y(r)}} \\{{t_{\tau}(r)} - {z(r)}}\end{bmatrix}}} & (15)\end{matrix}$

Referring to FIG. 2, surface S _(τ)(0) 209 is located at an altitude(t_(τ)=z_(τ) ^(O)) with respect to the vertical unit vector r _(t) _(τ)oriented upward. In the frame of the presentation of the volumedeformation and without loss of generality, z_(τ) ^(O) may be assumed tobe constant, e.g., equal to zero.

Referring to FIG. 2, surface S _(τ)(d) 207 is located at a distance (d)from S _(τ)(0) 209, implying that:

t _(τ)( r _(τ))=d−z _(τ) ^(O) ∀ r _(τ) ∈ S _(τ)(d)   (16)

such that:

$\begin{matrix}{{d < {{0{{\overset{\_}{S}}_{\tau}(d)}}\mspace{14mu} {is}\mspace{14mu} {located}\mspace{14mu} {above}\mspace{14mu} {{\overset{\_}{S}}_{\tau}(0)}}}{d > {{0{{\overset{\_}{S}}_{\tau}(d)}}\mspace{14mu} {is}\mspace{14mu} {located}\mspace{14mu} {above}\mspace{14mu} {{\overset{\_}{S}}_{\tau}(0)}}}} & (17)\end{matrix}$

FIG. 2 shows the folded present-day volume G_(τ) 202 resulting from thedeformation of restored volume G _(τ) 203 under tectonic forcesfollowing either a “minimal deformation” or a “flexural slip” tectonicstyle:

G _(τ)→tectonic forces→G_(τ)⊆G   (18)

Referring to FIG. 2, the following notation is used:

-   -   each horizontal plane S _(τ)(d) 207 is transformed into a curved        surface S_(τ)(d) 208 “parallel” (e.g., this notion of        “parallelism” may be characterized by equation (10)) to surface        H_(τ) 210=208 and each surface S_(τ)(d) 208 is a level set of        vertical restoration coordinate t_(τ)(r);    -   the images in G_(τ) 202 of the (u_(τ)), (v_(τ)) 205, 206 and        (t_(τ)) 204 coordinate axes initially rectilinear and contained        in volume G _(τ) 203 now consist of curvilinear coordinate axes        (223, 224) and 222.

As shown in FIG. 2, the part of the subsurface observed todaystratigraphically below H_(τ) 210 may be identified with the deformedvolume G_(τ) 202, e.g., such that:

-   -   S_(τ)(0) is assumed to be identical to the horizon H_(τ) 210 to        be restored:

S_(τ)(0)≡H_(τ)  (19)

-   -   which is equivalent to defining that, on horizon H_(τ) 216,        restored vertical coordinate t_(τ)(r) is equal to z_(τ) ^(O);    -   for any t<τ, the actual geologic horizon H_(τ) 216 is included        (212) into the deformed volume G_(τ) 202; note that, contrary to        surfaces {S_(τ)(d):d≥0} 208, horizons {H_(t):t<τ} 216 may be        non-parallel to {H_(τ)=S_(τ)(0)} 210=208;    -   after restoration of the volume G_(τ) 202 to its initial,        unfolded state G _(τ) 203:        -   all horizons, faults and geological objects included in            G_(τ) 202 are dragged up by the embedding volume            deformation,        -   H _(τ)≡S _(τ)(0)} 209=207 may be defined as the restored sea            floor as it was at geological time τ.

Minimizing Deformations and Volume Variation

With compaction handled separately in pre and post restoration steps,leaving aside the very particular case of clay and salt layers, tectonicforces generally induce no or negligible variations in volume.Therefore, restoration coordinates {u_(τ)(r), v_(τ)(r), t_(τ)(r)} may bechosen in such a way that the u_(τ)v_(τ)t_(τ)-transform 201 of thepresent-day volume G_(τ) 202 into the restored volume G _(τ) 203minimizes deformations and volume variations. This is achieved byconstraining restoration coordinates {u_(τ)(r), v_(τ)(r), t_(τ)(r)} tohonor the two following conditions in the present day G_(τ) domain:

-   -   Surfaces {S_(τ)(d):d≥0} 208 are level sets of the vertical        restoration coordinate t_(τ)(r) and, for any infinitely small        increment ε, the thickness of the thin slice of the volume        bounded by S_(τ)(d) and S_(τ)(d+ε) are, as much as possible,        constant and equal to ε. In other words, S_(τ)(d) and S_(τ)(d+ε)        are as parallel as possible. This is equivalent to honoring        equation (10) as precisely as possible.    -   In the frame of this invention, the consistency between the        depositional (e.g., GeoChron) model provided as input and its        restored version at geological time τ is preserved. Such a        consistency is preserved if, and only if, the uvt-transform and        the u_(τ)v_(τ)t_(τ)-transform of H_(τ) are identical. This is        achieved by honoring the following inventive boundary        conditions, referred to as the (u_(τ), v_(τ)) boundary        constraints:

$\begin{matrix}\left. {\forall{r_{\tau}^{o} \in H_{\tau}}} \middle| \begin{matrix}\left. 1 \right) & {{u_{\tau}\left( r_{\tau}^{o} \right)} = {u\left( r_{\tau}^{o} \right)}} \\\left. 2 \right) & {{v_{\tau}\left( r_{\tau}^{o} \right)} = {v\left( r_{\tau}^{o} \right)}}\end{matrix} \right. & (20) \\\left. {\forall{r_{\tau}^{o} \in H_{\tau}}} \middle| \begin{matrix}\left. 1 \right) & {{{grad}\mspace{14mu} {u_{\tau}\left( r_{\tau}^{o} \right)}} \simeq {{grad}\mspace{14mu} {u\left( r_{\tau}^{o} \right)}}} \\\left. 2 \right) & {{{grad}\mspace{14mu} {v_{\tau}\left( r_{\tau}^{o} \right)}} \simeq {{grad}\mspace{14mu} {v\left( r_{\tau}^{o} \right)}}}\end{matrix} \right. & (21)\end{matrix}$

-   -   Whilst taking the same given tectonic style into account        (minimal deformation or flexural slip) as the one honored by        paleo-geographic coordinates {u(r), v(r)}, lateral restoration        coordinates {u_(τ)(r), v_(τ)(r)} may be defined so that their        associated restoration deformations are minimized. To preserve        consistency with boundary conditions (20) and (21), this is        achieved by honoring the following inventive constraints:        -   in a minimal deformation tectonic style context, the            following “minimal deformation constraints” may be honored            by coordinates {u_(τ), v_(τ)}_(r) where t_(τ)(r) is given,            e.g., as follows:

$\begin{matrix}\left. {\forall{r \in G_{\tau}}} \middle| \begin{matrix}\left. 1 \right) & {\left\{ {{grad}\mspace{14mu} {u_{\tau} \cdot {grad}}\mspace{14mu} v_{\tau}} \right\}_{r} \simeq 0} \\\left. 2 \right) & {\left\{ {{grad}\mspace{14mu} {u_{\tau} \cdot {grad}}\mspace{14mu} v_{\tau}} \right\}_{r} \simeq 0} \\\left. 3 \right) & {\left\{ {{grad}\mspace{14mu} {u_{\tau} \cdot {grad}}\mspace{14mu} v_{\tau}} \right\}_{r} \simeq 0}\end{matrix} \right. & (22)\end{matrix}$

-   -   -   in a flexural slip tectonic style context, the following            “flexural slip constraint” is coupled (containing both            lateral restoration coordinates u_(τ) and v_(τ)) and may be            honored by coordinates {u_(τ), v_(τ)}_(r), e.g., as follows:

∀ r ∈ G_(τ): {grad_(S) u_(τ)·grad_(S) v_(τ)}_(r)≃0   (23)

-   -   -   where subscript “S” refers to a projection of the directions            of maximal change over iso-value surfaces of the restored            vertical coordinate tτ.

    -   So as not to conflict with equations (20) and (21), and contrary        to conventional depositional coordinates u and v (e.g., in the        GeoChron model), new constraints (22) and (23) do not constrain        ∥grad u_(τ)∥, ∥grad v_(τ)∥, ∥grad_(S) u_(τ)∥, or ∥grad_(S)        v_(τ)∥ to be equal to “1”.

Restoration

Referring to FIG. 1 and FIG. 2, at geological-time τ, the horizon H_(τ)210 to be restored was coincident with a given surface S_(τ)(0)(208=210) considered as the sea-floor. The task of restoration includes:

-   -   restoring horizon H_(τ) 210 to its initial, unfaulted and        unfolded state (e.g., mapping horizon H_(τ) onto the sea floor S        _(τ)(0)) 209 and    -   shifting all sedimentary terrains in such a way that, for each        point r ∈ G:        -   the particle of sediment currently located at point r moves            to its former, “restored” location, where the particle was            located at geological time τ,        -   no overlaps or voids/gaps are created in the subsurface.

At geological time τ, the sea floor S _(τ)(0) (209) is assumed to be acontinuous, unfaulted surface whose altitude z_(τ) ^(O) is a givenfunction z_(τ) ^(O) (u, v). In practice, S _(τ)(0) (209=207) istypically a flat, horizontal plane whose altitude z_(τ) ^(O) (u, v) atgeological time τ is constant. Accordingly, for concision, z_(τ) ^(O)may refer to a given function z_(τ) ^(O) (u(r), v(r)) which may or maynot be constant:

∀ r ∈ H_(τ): z_(τ) ^(O) stands for z_(r) ^(O)(u(r), v(r))   (24)

Compaction

Deformation of sedimentary terrains is typically induced both bytectonic forces and terrain compaction. In order to model separately theeffects of these phenomena, the restoration process may proceed asfollows:

-   -   First, in a preprocessing phase, a total decompaction may be        applied to the terrains to cancel the impact of compaction as it        is observed today, at the present day or current geological        time;    -   Next, the effects of compaction being canceled, a        depositional-based restoration process taking only tectonic        deformations into account (and not compaction) is applied to        restore the geometry of the subsurface as it would have been        observed at geological time τ;    -   Finally, in a post-processing phase, a partial recompaction is        applied to the restored terrains in order to take compaction        into account, as it could have been observed at geological time        τ.

Depositional Based Restoration

As an input to the restoration process, a given depositional (e.g.,GeoChron) model may be received from storage in a digital device (e.g.,from memory 150 of FIG. 15).

Referring to FIG. 2, a geological time τ may be selected that isassociated with the given horizon H_(τ) 210 to be restored and the givenaltitude z_(τ) ^(O) of the surface S _(τ) (0) 209 onto which the horizonH_(τ) should be restored.

The region G_(τ) 202 may be retrieved as the part of the depositionalmodel where geological time of deposition t(r) is less than or equal toτ (subsurface regions deposited in a layer deeper than or equal to thelayer deposited at time τ).

The set of faults may be split into a subset of τ-active faults cuttingH_(τ) 210 and a subset of τ-inactive faults which do not cut H_(τ).

A geologist or other user may decide to manually transfer some faultsfrom the τ-inactive fault set to the τ-active set, or vice versa, whichtypically causes greater restoration deformations. For example, manuallyaltering the set of automatically computed τ-active and τ-inactivefaults typically makes the restoration less accurate.

Four new 3D piecewise continuous discrete functions {u_(τ), v_(τ),t_(τ), ε_(τ)}_(r) may be created that are defined on grid Γ 100 whosediscontinuities occur only across τ-active faults;

Referring to FIG. 3, to remove discontinuities of discrete functions{u_(τ), v_(τ), t_(τ), ε_(τ)}_(r) across τ-inactive faults, for allτ-inactive faults F 300, one or more of the following inventive (e.g.,DSI) constraints may be installed on Γ 100, e.g., as:

$\begin{matrix}{\left. \begin{matrix}\left. 1 \right) & {{u_{\tau}\left( r_{F}^{\oplus} \right)} = {u_{\tau}\left( r_{F}^{\ominus} \right)}} \\\left. 2 \right) & {{v_{\tau}\left( r_{F}^{\oplus} \right)} = {v_{\tau}\left( r_{F}^{\ominus} \right)}} \\\left. 3 \right) & {{t_{\tau}\left( r_{F}^{\oplus} \right)} = {t_{\tau}\left( r_{F}^{\ominus} \right)}}\end{matrix} \right\} {\forall\left( {r_{F}^{\oplus},r_{F}^{\ominus}} \right)_{\tau}}} & (25) \\{\left. \begin{matrix}\left. 1 \right) & {{{grad}\mspace{14mu} {u_{\tau}\left( r_{F}^{\oplus} \right)}} = {{grad}\mspace{14mu} {u_{\tau}\left( r_{F}^{\ominus} \right)}}} \\\left. 2 \right) & {{{grad}\mspace{14mu} {v_{\tau}\left( r_{F}^{\oplus} \right)}} = {{grad}\mspace{14mu} {v_{\tau}\left( r_{F}^{\ominus} \right)}}} \\\left. 3 \right) & {{{grad}\mspace{14mu} {t_{\tau}\left( r_{F}^{\oplus} \right)}} = {{grad}\mspace{14mu} {t_{\tau}\left( r_{F}^{\ominus} \right)}}}\end{matrix} \right\} {\forall\left( {r_{F}^{\oplus},r_{F}^{\ominus}} \right)_{\tau}}} & (26) \\{\left. \begin{matrix}\left. 1 \right) & {{ɛ_{\tau}\left( r_{F}^{\oplus} \right)} = {ɛ_{\tau}\left( r_{F}^{\ominus} \right)}} \\\left. 2 \right) & {{{grad}\mspace{14mu} {ɛ_{\tau}\left( r_{F}^{\oplus} \right)}} = {{grad}\mspace{14mu} {ɛ_{\tau}\left( r_{F}^{\ominus} \right)}}}\end{matrix} \right\} {\forall\left( {r_{F}^{\oplus},r_{F}^{\ominus}} \right)_{\tau}}} & \begin{matrix}(27) \\(28)\end{matrix}\end{matrix}$

where (r_(F) ^(⊕),r_(F) ^(⊕))_(τ) (304,306) represents a pair of“mate-points” collocated on both sides of F 300 and assigned to F⁺ 103and F⁻ 104, respectively, and ε_(τ)(r) represents an error correctionconstraint. Constraints (25), (26), (27) and (28) may be referred tocollectively as “fault transparency constraints.”

Assuming that TH_(min)>0 is a given threshold chosen by a geologist orother user, fault transparency constraints (25), (26), (27) and (28) maybe locally installed at any point r_(F) on a τ-active fault F wherefault throw is lower than TH_(min). As a non-limitative example,TH_(min) may be equal to 1 meter.

Two new discrete vector fields r* and R_(τ) may be defined on 3D grid Γ100.

For each node α ∈ Γ 107:

-   -   r*(α) may be initialized as the initial location of α:

r*(α)=r(α)   (29)

-   -   a decompaction transformation C⁻(r) known in the art may be used        to move α vertically downward from its current compacted        altitude z(α) to a new decompacted (e.g., deeper) altitude:

r(α)←C ⁻(r(α))   (30)

Vertical Restoration Coordinate t_(τ)(r)

Referring to FIG. 2 and FIG. 7, the depositional uvt-transform 700 ofG_(τ) 202 is typically correct when equation (3) is honored. Based onthis observation, embodiments of the present invention adapt equation(3) for the inventive restoration technique, replacing the verticaldepositional coordinate t(r) by a vertical restoration coordinatet_(τ)(r) and replacing equation (3) by the following inventivethickness-preserving constraint to ensure layer thickness is preservedand surfaces {S_(τ)(d):d≥0} are parallel:

∥grad t_(τ)(r)∥≃1 ∀ r ∈ G_(τ)  (31)

In addition, to allow H_(τ) 210 to be restored on surface S _(τ)(0) 209,the vertical restoration coordinate t_(τ)(r) may honor the followingboundary condition, e.g., as a DSI constraint on grid Γ 100, referred toas the “paleo-sea-floor constraint”:

t _(τ)(r _(H))=z _(τ) ^(O) ∀ r _(H) ∈ H _(τ)  (32)

Due to its non-linearity, thickness-preserving equation (31) cannot beimplemented as a DSI constraint, which must be linear. In order toincorporate the thickness-preserving equation into the restoration modelusing the DSI method, various linear surrogates of equation (31) may beused to approximate t_(τ)(r) as follows:

-   -   Referring to FIG. 1, to approximate thickness-preserving        equation (31), as a non-limitative example, install the        following DSI constraints on the grid Γ 100 where r_(T⋄) and        r_(T)* are arbitrary points belonging to any pair (T^(⋄), T*) of        adjacent cells 108 of grid Γ 100 and where N(r_(h)) is the field        of unit vectors defined on H_(τ), orthogonal to H_(τ) and        oriented in the direction of younger terrains:

grad t _(τ)(r _(H))=N(r _(H)) • r _(H) ∈ H _(τ)  1)

grad t_(τ)(r_(T⋄))≃grad t_(τ)(r_(T)*) ∀ (T^(⋄), T*)   2)

(33)

-   -   Referring to FIG. 1, to approximate thickness-preserving        equation (31), alternatively as a non-limitative example,        install the inventive DSI constraints on grid Γ 100 as follows:

  (34)

$\begin{matrix}\left. 1 \right) & {{{grad}\mspace{14mu} {t_{\tau}(r)}} = \frac{{grad}\mspace{14mu} {t(r)}}{{{grad}\mspace{14mu} {t(r)}}}} & {\forall{r \in G_{\tau}}} \\\left. 2 \right) & {{{grad}\mspace{14mu} {t_{\tau}\left( r_{T\; ♦} \right)}} \simeq {{grad}\mspace{14mu} {t_{\tau}\left( r_{T^{*}} \right)}}} & {\forall\left( {T^{♦},T^{*}} \right)}\end{matrix}$

where r_(T⋄) and r_(T)* are arbitrary points belonging to any pair(T^(⋄), T*) of adjacent cells of grid Γ 100 (e.g., the centers of T^(⋄)and T*, respectively).

Constraints (33) and (34) are only examples of possiblesurrogate-thickness-preserving constraints. Other examples of suchsurrogate thickness-preserving constraints may be used.

Referring to FIG. 3, contrary to constraint (33)-(1), new inventiveconstraint (34)-(1) benefits from the geologic observation that,throughout the entire domain G_(τ) 202, surfaces {S_(τ)(d):d≥0} 208generally have a shape roughly similar to the level sets of the geologictime of deposition t(r);

Assuming that constraints (32) and (33) or (34) are installed on grid Γ100, a first approximation of vertical restoration coordinate t′_(τ)(r)may be computed by running the DSI method on grid Γ 100.

Honoring constraint (31) significantly increases the accuracy of therestoration model and a violation of this constraint not only degradesthe accuracy of the vertical restoration coordinate t_(τ)(r) but alsothe horizontal restoration coordinates {u_(τ)(r), v_(τ)(r)} as they arelinked to t_(τ)(r) (e.g., by equations (22) and (23)). Accordingly,there is a great need for validating any approximation technique used tocompute t_(τ)(r).

To test the accuracy of the various approximations of t_(τ)(r), anexample geological terrain is provided in FIG. 4. Despite the apparentsimplicity of this terrain, because the thicknesses of the layers vary,this test example is challenging and useful in comparing the accuracy ofinventive embodiments with other conventional techniques.

FIG. 5 shows histograms 501 and 502 of the distributions of ∥gradt_(τ)∥, where t_(τ) is approximated using constraints (33) or (34),respectively, in the example geological terrain G_(τ) 202 shown in FIG.4. FIG. 5 shows that when t_(τ) is approximated by constraints (33) or(34), ∥grad t_(τ)∥ significantly differs from “1”. Therefore, whileconstraints (34) provide a better approximation of thethickness-preserving equation (31) than constraints (33), both of theseapproximations are inaccurate.

Similarly, FIG. 9 shows histograms 901 and 902 of relative variations ofvolume ΔV/V induced by the restoration of H_(τ) 210 over G_(τ) 202 shownin FIG. 4, where t_(τ) is approximated using constraints (33) or (34),respectively. Ideally, a restoration transformation should minimizevariations in volume ΔV/V from the present day to the restored model.FIG. 9 however shows that a restoration based on constraints (33) or(34) results in a volume variation ΔV/V that significantly differs fromthe ideal value of “0”. While constraints (34) result in a smallervolume variation ΔV/V than constraints (33), both of theseapproximations induce a significant volume variation ΔV/V and areinaccurate.

Improving Vertical Restoration Coordinate t_(τ)(r)

An approximation of the vertical restoration coordinate t′_(τ)(r) may beimproved by a “t_(τ)-incremental improvement” constraint, e.g., asfollows:

t _(τ)(r)=t′ _(r)(r)+ε_(τ)(r) ∀ r ∈ G _(τ)  (35)

where ε_(τ)(r) is an error correction term, e.g., as characterizedbelow.

Accordingly, assuming that an initial approximation t′_(τ)(r) hasalready been obtained, to compute an improved version of t_(τ)(r), thefollowing inventive incremental procedure may be executed:

-   -   For each point r_(H) ∈ H_(τ), set the following equation as an        inventive sea-floor-error constraint e.g., using the DSI method:

ε_(τ)(r _(H))=0 ∀ r _(H) ∈ H _(τ)  (36)

-   -   this constraint assumes that constraint (32) remains honored.    -   For each cell T ∈ Γ ∩ G_(τ) 108, choose a point r_(τ) in the        cell T (e.g., its center) and install the new linear        thickness-preserving constraint, e.g., using the DSI method as        follows:

grad ε_(τ)(r_(τ))·grad t′_(τ)(r_(τ))≃½{1−∥grad t′τ(r_(τ))∥²}  (37)

-   -   This constraint is linear, deduced from a linear second order        approximation of equation (31). Further, this constraint ensures        that, after applying the t_(τ)-incremental improvement        correction constraint (35), the local value of ∥grad t_(τ)(r)∥        at any point r ∈ G_(τ) is as close as possible to “1.”    -   For each sampling point r located on a τ-active fault, install        for ε_(τ)(r) the following inventive DSI constraint referred to        as the “t_(τ)-incremental boundary” constraint:

grad ε_(τ)(r)×grad t_(τ)(r)≃0   (38)

-   -   This constraint specifies that, after applying correction        constraint (35), in the close neighborhood of τ-active faults,        the shape of level sets of t_(τ)(r) remains roughly unchanged.    -   To ensure piecewise continuity of the error correction ε_(τ)(r),        install DSI gradient smoothness constraints, known in the art,        for the error correction ε_(τ)(r).    -   Assuming that constraints (36), (37) and (38) are installed on        grid Γ 100, to interpolate the error correction ε_(τ)(r), run        DSI on grid Γ 100.    -   For each node α ∈ Γ 107, update the vertical restoration        coordinate t_(τ)(α) as follows:

t _(τ)(α)=t′ _(r)(α)+ε_(τ)(α)   (39)

-   -   In the test case represented by FIG. 4 and FIG. 5, the histogram        503 of the distribution of ∥grad t_(τ)∥, where t_(τ) is        approximated by constraints (37) over G_(τ) 202 is now        considerably better than histograms 501 and 502 obtained with        constraints (33) or (34), respectively. In particular:        -   As specified by equation (10), distribution 503 is now            centered on value “1”. This condition is of paramount            importance to minimize deformations during the restoration            process generated by a u_(τ)v_(τ)t_(τ)-transform.        -   The standard deviation of distribution 503 is considerably            reduced as compared to the relatively wider standard            deviation of distributions 501 and 502.    -   Moreover, in the test case represented in FIG. 4 and FIG. 9, the        histogram 903 of the distribution of relative variations of        volume ΔV/V induced by a restoration of H_(τ) 210 over G_(τ)        202, where t_(τ) is obtained using constraints (37) is        significantly more accurate than histograms 901 and 902, where        t_(τ) is obtained using constraints (33) or (34), respectively.        The center of histogram 903 of volume variation ΔV/V is closer        to the ideal value “0” than histograms 901 and 902, which        indicates that variations of volume are better minimized after        applying second order constraints (37) than constraints (33) or        (34).

Horizontal Restoration Coordinates {u_(τ)(r), v_(τ)(r)}

Referring to FIG. 2, with respect to surfaces {S_(τ)(d):d≥0} 208,horizontal restoration coordinates {u_(τ)(r), v_(τ)(r)} play a rolesimilar to the one played by paleo-geographic coordinates {u(r), v(r)}with respect to horizons {H_(t):t≥0} 216 of the depositional modelprovided as input. Based on this similarity, horizontal restorationcoordinates {u_(τ)(r), v_(τ)(r)} may be generated as follows:

-   -   install equations (20) and (21) as inventive boundary        constraints.    -   for all points r ∈ G_(τ) 214, define as follows inventive        vectors fields a_(τ)(r) and b_(τ)(r) respectively, referred to        as the “τ-axe” and “τ-coaxe” vector fields:

a _(τ)(r)=grad t _(τ)(r)×grad u(r)×grad t _(τ)(r)

b _(τ)(r)=grad t _(τ)(r)×a _(τ)(r)   (40)

The τ-axe and τ-coaxe vector fields a_(τ)(r) and b_(τ)(r) differconsiderably from the local axe and co-axe vectors fields a(r) and b(r),e.g., as discussed in U.S. Pat. No. 8,711,140, which is incorporated byreference herein in its entirety. These new τ-axe and τ-coaxe vectorsa_(τ)(r) and b_(τ)(r) strongly depend on the new vertical restorationcoordinate t_(τ)(r) (e.g., already computed as above) and also take intoaccount the gradient of the paleo-geographic coordinate u(r) (e.g.,associated to the depositional model provided as input).

-   -   if the tectonic style is minimal deformation then, to        approximate equations (22), install the following inventive        “surrogate minimal-deformation” constraints e.g., using the DSI        method:

$\begin{matrix}{\forall{r \in {G_{\tau}:\left| \begin{matrix}{{{grad}\mspace{14mu} {{u_{\tau}(r)} \cdot {b_{\tau}(r)}}} \simeq 0} \\{{{grad}\mspace{14mu} {{v_{\tau}(r)} \cdot {\alpha_{\tau}(r)}}} \simeq 0}\end{matrix} \right.}}} & (41)\end{matrix}$

-   -   if the tectonic style is flexural slip then, to approximate        equations (23), install the following inventive “surrogate        flexural-slip” constraints e.g., using the DSI method:

$\begin{matrix}{\forall{r \in {G_{\tau}:\left| \begin{matrix}{{{grad}_{s}\mspace{14mu} {{u_{\tau}(r)} \cdot {b_{\tau}(r)}}} \simeq 0} \\{{{grad}_{s}\mspace{14mu} {{v_{\tau}(r)} \cdot {\alpha_{\tau}(r)}}} \simeq 0}\end{matrix} \right.}}} & (42)\end{matrix}$

-   -   -   where subscript “S” refers to a projection of the directions            of maximal change over iso-value surfaces of the restored            vertical coordinate tτ.

    -   Referring to FIG. 1 and FIG. 6, to prevent the        u_(τ)v_(τ)t_(τ)-transform used as a restoration operator from        generating gaps and overlaps along τ-active faults, specific        constraints may be added along fault striae induced by        twin-points of the depositional model provided as input. For        that purpose, for each pair of twin-points (r_(F) ⁺,r_(F) ⁻)        (101, 102) located on faces F⁺ 103 and F⁻ 104 of a τ-active        fault F 105, respectively, a process may proceed according to as        follows:        -   retrieve the fault stria σ(r_(F) ⁻) 600 passing through twin            points (r_(F) ⁺, r_(F) ⁻) (101, 102), and        -   on curve σ(r_(F) ⁻) 600, retrieve a point {tilde over            (r)}_(F) ⁺ 601 located on F⁺ 103 and such that t_(τ)({tilde            over (r)}_(F) ⁺) is approximately equal to t_(τ)(r_(F) ⁻).

    -   install the following inventive “τ-fault-striae” constraints        e.g., using the DSI method:

$\begin{matrix}\left| \begin{matrix}\left. 1 \right) & {{u_{\tau}\left( {\overset{\sim}{r}}_{F}^{+} \right)} \simeq {u_{\tau}\left( r_{F}^{-} \right)}} \\\left. 2 \right) & {{v_{\tau}\left( {\overset{\sim}{r}}_{F}^{+} \right)} \simeq {v_{\tau}\left( r_{F}^{-} \right)}}\end{matrix} \right. & (43)\end{matrix}$

-   -   To ensure piecewise continuity of horizontal restoration        coordinates u_(τ)(r) and v_(τ)(r), install gradient smoothness        constraints e.g., using the DSI method.    -   To compute the pair of horizontal restoration coordinates        {u_(τ)(r), v_(τ)(r)} honoring constraints (20), (41 or 42) and        (43), run DSI on grid Γ 100.

Computing the Restoration R_(τ)(r)

The restoration vector field R_(τ)(r) represents the field ofdeformation vectors from the present day (e.g., xyz) space to therestoration (e.g., u_(τ)v_(τ)t_(τ)) space, e.g., computed from theu_(τ)v_(τ)t_(τ)-transform.

Referring to FIG. 1, for each node α 107 of 3D grid Γ 100, move a torestored location r(α), e.g., defined as follows:

r (α)=u _(τ)(α)·r _(x) +v _(τ)(α)·r _(y) +t _(τ)(α)·r _(z)   (44)

For each node α 107 of 3D grid Γ 100:

-   -   if, to compute vertical restoration coordinate t_(τ)(r),        compaction was taken into account, then, using a recompaction        operator C⁺(r) known in the art, move α vertically upward from        its current decompacted altitude z(α) to a new recompacted        (shallower) altitude:

r(α)←C ⁺(r(α))   (45)

-   -   save the restoration vector Rτ(α) on a digital device:

Rτ(α)=r(α)·r*(α)   (46)

-   -   where r*(α) is defined e.g., in equation (29).    -   reset location r(α) of α to its initial location before        restoration:

r(α)←r*(α)   (47)

-   -   stop.

Scanning the Subsurface Through Time

Consider a series of geological restoration times {τ₁<τ₂< . . . <τ_(n)}associated with reference horizons H_(τ) ₁ , H_(τ) ₂ , . . . , H_(τ)_(n) , respectively. Using the restoration method described herein, foreach (τ_(i)=τ), build and store on a digital device a restoration vectorfield R_(τ) _(i) (r)=R_(τ)(r), e.g., as:

$\begin{matrix}{\mspace{79mu} {\begin{matrix}\tau_{1} & < & \tau_{2} & < & \cdots & < & \tau_{n} \\ \updownarrow & \; & \updownarrow & \; & \; & \; & \updownarrow \\R_{\tau_{\text{?}}} & \; & R_{\tau_{2}} & \; & \; & \; & R_{\tau_{n}}\end{matrix}{\text{?}\text{indicates text missing or illegible when filed}}}} & (48)\end{matrix}$

In addition to these reference restoration times, an additionalrestoration time τ_(n+1) may be added to be associated with thehorizontal plane H_(t) _(n+1) located at a constant altitude z_(τ)_(n+1) ⁰ of the sea level. Time τ_(n+1) may be the present daygeological time and, provided that τ_(n+1) is greater than τ_(n), anyarbitrary value may be chosen for τ_(n+1). As a non-limitative example,τ_(n+1) may be defined as:

τ_(n+1)=τ_(n)+1   (49)

Because τ_(n+1) is the present day, applying the restoration vectorfield Rτ_(n+1)(r) to the present day horizon H_(t) _(n+1) should leaveH_(t) _(n+1) unchanged e.g., as follows:

Rτ _(n+1)(r)=0 ∀ r ∈ G   (50)

To explore subsurface evolution throughout geological times, a processmay proceed as follows:

-   -   as input, read a depositional (e.g., GeoChron) model and        associated series of restoration vector fields {R_(τ) ₁ , R_(τ)        ₂ , . . . , R_(τ) _(n+1) } stored on a digital device;    -   using an input device such as, in a non-limitative example, the        keyboard of a computer or the wheel of a computer mouse, select        a geological time τ_(i) in the given list of geological times        {τ₁<τ₂< . . . <τ_(n+1)};    -   for each vertex α ∈ Γ 107, save a copy r*(α) of the location of        this node in the depositional model given as input;    -   apply the restoration vector field R_(τ) _(i) (r) to the        depositional model given as input;    -   display the restored model on a device such as, in a        non-limitative example, a display (e.g., display 180 of FIG.        15), such as, a screen, a hologram or a 3D printer;    -   repeat the previous operations as long as desired.    -   optionally, to modify the geometry of the horizons at geological        time τ_(i), use a computerized tool known in the art to edit the        geological time of deposition t(r);    -   for each vertex α 107 of 3D grid Γ 100, use copy r*(α) to        restore r(α) to its present day location:

r(α)←r*(α) ∀ α ∈ Γ  (51)

-   -   such an operation implicitly and automatically propagates the        modifications of the geometry of horizons optionally performed        above;    -   return to the first step above.

Geological Tomography

Geological models are generated using geological or seismic tomographytechnology. Geological tomography generates an image of the interiorsubsurface of the Earth based on geological data collected bytransmitting a series of incident waves and receiving reflections ofthose waves across discontinuities in the subsurface. A transmitter maytransmit signals, for example, acoustic waves, compression waves orother energy rays or waves, that may travel through subsurfacestructures. The transmitted signals may become incident signals that areincident to subsurface structures. The incident signals may reflect atvarious transition zones or geological discontinuities throughout thesubsurface structures, such as, faults or horizons. The reflectedsignals may include seismic events. A receiver may collect data, forexample, reflected seismic events. The data may be sent to a modelingmechanism that may include, for example, a data processing mechanism andan imaging mechanism.

Reference is made to FIG. 14, which is a schematic illustration of ageological tomography technique in which a series of incident rays 111and reflected rays 121 are propagated through a subsurface region of theEarth 30 to image the subsurface, according to an embodiment of theinvention.

One or more transmitter(s) (e.g., 190 of FIG. 15) located at incidentlocation(s) 60 may emit a series of incident rays 111. Incident rays 111may include for example a plurality of energy rays related to signalwaves, e.g., sonic waves, seismic waves, compression waves, etc.Incident rays 111 may be incident on, and reflect off of, a subsurfacestructure or surface 90 at a reflection point 50. Multiple reflectionpoints 50 may be identified or imaged or displayed in conjunction todisplay, for example, a horizon.

One or more receiver(s) (e.g., 140 of FIG. 15) located at reflectedlocation(s) 65 may receive the reflection rays 121. Reflection rays 121may be the reflected images of incident rays 111, for example, afterreflecting off of image surface 90 at target point 50. The angle ofreflection 55 may be the angle between corresponding incident rays 111and reflected rays 121 at reflection point 50. An incident rays 111 anda corresponding reflected rays 121 may propagate through a cross-sectionof a subsurface structure 30. Incident rays 111 may reflect off of asubsurface feature 90 at a reflection point 50, for example, a point onan underground horizon, the seafloor, an underground aquifer, etc.

One or more processor(s) (e.g., 140 of FIG. 15) may reconstituteincident and reflected rays 111 and 121 to generate an image thesubsurface 30 using an imaging mechanism. For example, a commonreflection angle migration (CRAM) imaging mechanism may image reflectionpoints 50 by aggregating all reflected signals that may correspond to areflection point, for example, reflected signals that may have the samereflection angle. In other examples, imaging mechanisms may aggregatereflected signals that may have the same reflection offset (distancebetween transmitter and receiver), travel time, or other suitableconditions.

The processor(s) may compose all of the reflection points 50 to generatean image or model of the present day underground subsurface of the Earth30. The processor(s) may execute a restoration transformation (e.g.,u_(τ)v_(τ)t_(τ)-transform 201) to transform the present day model ofsubsurface 30 to a restored subsurface image 203 at a restorationtime_(τ). One or more display(s) (e.g., 180 of FIG. 15) may visualizethe present day subsurface image 30 and/or the restored subsurface image203.

System Overview

Reference is made to FIG. 15, which schematically illustrates a systemincluding one or more transmitter(s), one or more receiver(s) and acomputing system in accordance with an embodiment of the presentinvention. Methods disclosed herein may be performed using a system 1505of FIG. 15.

System 1505 may include one or more transmitter(s) 190, one or morereceiver(s) 120, a computing system 130, and a display 180. Theaforementioned data, e.g., seismic data used to form intermediate dataand finally to model subsurface regions, may be ascertained byprocessing data generated by transmitter 190 and received by receiver120. Intermediate data may be stored in memory 150 or other storageunits. The aforementioned processes described herein may be performed bysoftware 160 being executed by processor 140 manipulating the data.

Transmitter 190 may transmit signals, for example, acoustic waves,compression waves or other energy rays or waves, that may travel throughsubsurface (e.g., below land or sea level) structures. The transmittedsignals may become incident signals that are incident to subsurfacestructures. The incident signals may reflect at various transition zonesor geological discontinuities throughout the subsurface structures. Thereflected signals may include seismic data.

Receiver 120 may accept reflected signal(s) that correspond or relate toincident signals, sent by transmitter 190. Transmitter 190 may transmitoutput signals. The output of the seismic signals by transmitter 190 maybe controlled by a computing system, e.g., computing system 130 oranother computing system separate from or internal to transmitter 190.An instruction or command in a computing system may cause transmitter190 to transmit output signals. The instruction may include directionsfor signal properties of the transmitted output signals (e.g., such aswavelength and intensity). The instruction to control the output of theseismic signals may be programmed in an external device or program, forexample, a computing system, or into transmitter 190 itself.

Computing system 130 may include, for example, any suitable processingsystem, computing system, computing device, processing device, computer,processor, or the like, and may be implemented using any suitablecombination of hardware and/or software. Computing system 130 mayinclude for example one or more processor(s) 140, memory 150 andsoftware 160. Data 155 generated by reflected signals, received byreceiver 120, may be transferred, for example, to computing system 130.The data may be stored in the receiver 120 as for example digitalinformation and transferred to computing system 130 by uploading,copying or transmitting the digital information. Processor 140 maycommunicate with computing system 130 via wired or wireless command andexecution signals.

Memory 150 may include cache memory, long term memory such as a harddrive, and/or external memory, for example, including random accessmemory (RAM), read only memory (ROM), dynamic RAM (DRAM), synchronousDRAM (SD-RAM), flash memory, volatile memory, non-volatile memory, cachememory, buffer, short term memory unit, long term memory unit, or othersuitable memory units or storage units. Memory 150 may storeinstructions (e.g., software 160) and data 155 to execute embodiments ofthe aforementioned methods, steps and functionality (e.g., in long termmemory, such as a hard drive). Data 155 may include, for example, rawseismic data collected by receiver 120, instructions for building a mesh(e.g., 100), instructions for partitioning a mesh, and instructions forprocessing the collected data to generate a model, or other instructionsor data. Memory 150 may also store instructions to divide and modelτ-active faults and τ-inactive faults. Memory 150 may generate and storethe aforementioned constraints, restoration transformation (e.g.,u_(τ)v_(τ)t_(τ)-transform 201), restoration coordinates (e.g., u_(τ),v_(τ), t_(τ)), a geological-time and paleo-geographic coordinates (e.g.,u, v, t), a model representing a structure when it was originallydeposited (e.g., in uvt-space), a model representing a structure at anintermediate restoration time (e.g., in u_(τ), v_(τ), t_(τ)-space),and/or a model representing the corresponding present day structure in acurrent time period (e.g., in xyz-space). Memory 150 may store cells,nodes, voxels, etc., associated with the model and the model mesh.Memory 150 may also store forward and/or reverse u_(τ), v_(τ),t_(τ)-transformations to restore present day models (e.g., in xyz-space)to restored models (e.g., in u_(τ), v_(τ), t_(τ)-space), and vice versa.Memory 150 may also store the three-dimensional restoration vectorfields, which when applied to the nodes of the initial present daymodel, move the nodes of the initial model to generate one of theplurality of restored models. Applying a restoration vector field tocorresponding nodes of the present day model may cause the nodes to“move”, “slide”, or “rotate”, thereby transforming modeled geologicalfeatures represented by nodes and cells of the initial model. Data 155may also include intermediate data generated by these processes and datato be visualized, such as data representing graphical models to bedisplayed to a user. Memory 150 may store intermediate data. System 130may include cache memory which may include data duplicating originalvalues stored elsewhere or computed earlier, where the original data maybe relatively more expensive to fetch (e.g., due to longer access time)or to compute, compared to the cost of reading the cache memory. Cachememory may include pages, memory lines, or other suitable structures.Additional or other suitable memory may be used.

Computing system 130 may include a computing module havingmachine-executable instructions. The instructions may include, forexample, a data processing mechanism (including, for example,embodiments of methods described herein) and a modeling mechanism. Theseinstructions may be used to cause processor 140 using associatedsoftware 160 modules programmed with the instructions to perform theoperations described. Alternatively, the operations may be performed byspecific hardware that may contain hardwired logic for performing theoperations, or by any combination of programmed computer components andcustom hardware components.

Embodiments of the invention may include an article such as anon-transitory computer or processor readable medium, or a computer orprocessor storage medium, such as for example a memory, a disk drive, ora USB flash memory, encoding, including or storing instructions, e.g.,computer-executable instructions, which when executed by a processor orcontroller, carry out methods disclosed herein.

Display 180 may display data from transmitter 190, receiver 120, orcomputing system 130 or any other suitable systems, devices, orprograms, for example, an imaging program or a transmitter or receivertracking device. Display 180 may include one or more inputs or outputsfor displaying data from multiple data sources or to multiple displays.For example, display 180 may display visualizations of subsurface modelsincluding subsurface features, such as faults, horizons andunconformities, as a present day subsurface image (e.g., 202), arestored subsurface image (e.g., 203) and/or a depositional model (e.g.,703). Display 180 may display one or more present day model(s),depositional model(s), restoration model(s), as well as a series ofchronologically sequential restoration models associated with a sequenceof respective restoration times (e.g., τ₁<τ₂<τ₃<τ₄, as shown in FIG.13). The models may be displayed one at a time, two at a time, or manyat a time (e.g., the number selected by a user or automatically based onthe difference between models or the total number of models). Display180 may display the models in a sequence of adjacent models, throughwhich a user may scan (e.g., by clicking a ‘next’ or ‘previous’ buttonwith a pointing device such as a mouse or by scrolling through themodels).

Input device(s) 165 may include a keyboard, pointing device (e.g.,mouse, trackball, pen, touch screen), or cursor direction keys, forcommunicating information and command selections to processor 140. Inputdevice 165 may communicate user direction information and commandselections to the processor 140. For example, a user may use inputdevice 165 to select one or more preferred models from among theplurality of perturbed models, recategorize faults as τ-active faultsand τ-inactive, or edit, add or delete subsurface structures.

Processor 140 may include, for example, one or more processors,controllers or central processing units (“CPUs”). Software 160 may bestored, for example, in memory 150. Software 160 may include anysuitable software, for example, DSI software.

Processor 140 may generate a present day subsurface image (e.g., 202), arestored subsurface image (e.g., 203) and/or a depositional model (e.g.,703), for example, using data 155 from memory 150. In one embodiment, amodel may simulate structural, spatial or geological properties of asubsurface region, such as, porosity or permeability through geologicalterrains.

Processor 140 may initially generate a three dimensional mesh, lattice,grid or collection of nodes (e.g., 100) that spans or covers a domain ofinterest. The domain may cover a portion or entirety of thethree-dimensional subsurface region being modeled. Processor 140 mayautomatically compute the domain to be modeled and the correspondingmesh based on the collected seismic data so that the mesh covers aportion or the entirety of the three-dimensional subsurface region fromwhich geological data is collected (e.g., the studied subsurfaceregion). Alternatively or additionally, the domain or mesh may beselected or modified by a user, for example, entering coordinates orhighlighting regions of a simulated optional domain or mesh. Forexample, the user may select a domain or mesh to model a region of theEarth that is greater than a user-selected subsurface distance (e.g.,100 meters) below the Earth's surface, a domain that occurs relative togeological features (e.g., to one side of a known fault or riverbed), ora domain that occurs relative to modeled structures (e.g., betweenmodeled horizons H(t₁) and H(t₁₀₀)). Processor 140 may execute software160 to partition the mesh or domain into a plurality ofthree-dimensional (3D) cells, columns, or other modeled data (e.g.,represented by voxels, pixels, data points, bits and bytes, computercode or functions stored in memory 150). The cells or voxels may havehexahedral, tetrahedral, or any other polygonal shapes, and preferablythree-dimensional shapes. Alternatively, data may includezero-dimensional nodes, one-dimensional segments, two-dimensional facetand three-dimensional elements of volume, staggered in athree-dimensional space to form three-dimensional data structures, suchas cells, columns or voxels. The cells preferably conform to andapproximate the orientation of faults and unconformities. Each cell mayinclude faces, edges and/or vertices. Each cell or node may correspondto one or more particles of sediment in the Earth (e.g., a node mayinclude many cubic meters of earth, and thus many particles).

Data collected by receiver 120 after the time of deposition in a currentor present time period, include faults and unconformities that havedeveloped since the original time of deposition, e.g., based on tectonicmotion, erosion, or other environmental factors, may disrupt the regularstructure of the geological domain. Accordingly, an irregular mesh maybe used to model current geological structures, for example, so that atleast some faces, edges, or surfaces of cells are oriented parallel tofaults and unconformities, and are not intersected thereby. In oneembodiment, a mesh may be generated based on data collected by receiver120, alternatively, a generic mesh may be generated to span the domainand the data collected by receiver 120 may be used to modify thestructure thereof. For example, the data collected may be used togenerate a set of point values at “sampling point”. The values at thesepoints may reorient the nodes or cells of the mesh to generate a modelthat spatially or otherwise represents the geological data collectedfrom the Earth. Other or different structures, data points, or sequencesof steps may be used to process collected geological data to generate amodel. The various processes described herein (e.g., restoring ageological model using τ-active and τ-inactive faults, or restoring ageological model using a new thickness-preserving constraint) may beperformed by manipulating such modeling data.

Restoration coordinates may be defined at a finite number of nodes orsampling points based on real data corresponding to a subsurfacestructure, e.g., one or more particles or a volume of particles ofEarth. Restoration coordinates may be approximated between nodes tocontinuously represent the subsurface structure, or alternatively,depending on the resolution in which the data is modeled may representdiscrete or periodic subsurface structures, e.g., particles or volumesof Earth that are spaced from each other.

The computing system of FIG. 15 may accept the data used in theoperations of FIGS. 16, 17 and 19 as for example a set of data generatedby tomographic scanning of a subsurface geological region of the Earthas disclosed in reference to FIG. 14, or such data augmented by anotherprocess. The computing system may accept one or more of seismic and welldata. The computing device may generate one or more of seismic and welldata.

“Restoration” or “intermediate” time τ may refer to a time in the pastbefore the present day and after a time when an oldest or deepesthorizon surface in the 3D model was deposited. “Restoration” or“intermediate” transformation or model may refer to a model or image ofthe surface as it was configured at the “intermediate” time in the pastτ. An intermediate horizon may refer to a horizon that was deposited atthe “intermediate” time τ, which is located above the deepest horizonand below the shallowest horizon.

“Time” including the present-day, current or present time, the pastrestoration time τ, and/or the depositional time t, may refer togeological time periods that span a duration of time, such as, periodsof thousands or millions of years.

“Geological-time” t(r) may refer to the time of deposition when aparticle of sediment represented by point r was originally deposited inthe Earth. In some embodiments, the geological-time of the depositionmay be replaced, e.g., by any arbitrary monotonic increasing function ofthe actual geological-time. It is a convention to use an monotonicallyincreasing function, but similarly an arbitrary monotonic decreasingfunction may be used. The monotonic function may be referred to as the“pseudo-geological-time”.

The geological-time of the deposition and restoration time of particlesare predicted approximate positions since past configurations can nottypically be verified.

“Current” or “present day” location for a particle (or data structurerepresenting one or more particles) or subsurface feature may refer tothe location of the item in the present time, as it is measured.

In stratified terrain, layers, horizons, faults and unconformities maybe curvilinear surfaces which may be for example characterized asfollows.

-   -   A horizon, Hτ, may be a surface corresponding to a plurality of        particles of sediment which were deposited approximately at        substantially the same geological-time, τ.    -   A fault may be a surface of discontinuity of the horizons that        may have been induced by a relative displacement of terrains on        both sides of such surfaces. In other words, the geological-time        of deposition of the sediments is discontinuous across each        fault. Faults may cut horizons and may also cut other faults.    -   An unconformity may be a surface of discontinuity of the        horizons that may have been induced by for example an erosion of        old terrains replaced by new ones. In other words, similarly to        faults, the geological-time of deposition of the sediments is        discontinuous across each unconformity.

Terrain deformed in the neighborhood of a point r in the G-space mayoccur according to a “minimal deformation” tectonic style when, in thisneighborhood:

the strain tensor is approximately equal to the null tensor.

Terrain deformed in the neighborhood of a point r in the G-space mayoccur according to a “flexural slip” tectonic style when, in thisneighborhood:

-   -   the length of any small increment of distance d(r) on the        horizon passing through point r is preserved, e.g., in any        direction, and,    -   the volume of the terrains in the neighborhood of point r does        not vary.

Discrete-Smooth-Interpolation (DSI) is a method for interpolating orapproximating values of a function f(x,y,z) at nodes of a 3D grid ormesh Γ (e.g., 100), while honoring a given set of constraints. The DSImethod allows properties of structures to be modeled by embedding dataassociated therewith in a (e.g., 3D Euclidean) modeled space. Thefunction f(x,y,z) may be defined by values at the nodes of the mesh, Γ.The DSI method allows the values of f(x,y,z) to be computed at the nodesof the mesh, Γ, so that a set of one or more (e.g., linear) constraintsare satisfied. DSI generally only applies linear constraints on themodel.

In some embodiments, bold symbols represent vectors or multi-dimensional(e.g., 3D) functions or data structures.

In some embodiments, the “simeq” symbol “≃” or “≅” may meanapproximately equal to, e.g., within 1-10% of, or in a least squaressense. In some embodiments, the symbol “≡” may mean identical to, ordefined to be equal to.

While embodiments of the invention describe the input depositional modelas the GeoChron model, any other depositional model visualizing thepredicted configuration of each particle, region or layer at itsrespective time of depositional may be used.

While embodiments of the invention describe the present day coordinatesas xyz, the restoration coordinates as u_(τ)v_(τ)t_(τ), the depositionalcoordinates as uvt, the restoration transformation as au_(τ)v_(τ)t_(τ)-transform, and the depositional transformation as auvt-transform, any other coordinates or transformations may be used.

In the foregoing description, various aspects of the present inventionhave been described. For purposes of explanation, specificconfigurations and details have been set forth in order to provide athorough understanding of the present invention. However, it will alsobe apparent to one skilled in the art that the present invention may bepracticed without the specific details presented herein. Furthermore,well known features may have been omitted or simplified in order not toobscure the present invention. Unless specifically stated otherwise, asapparent from the following discussions, it is appreciated thatthroughout the specification discussions utilizing terms such as“processing,” “computing,” “calculating,” “determining,” or the like,refer to the action and/or processes of a computer or computing system,or similar electronic computing device, that manipulates and/ortransforms data represented as physical, such as electronic, quantitieswithin the computing system's registers and/or memories into other datasimilarly represented as physical quantities within the computingsystem's memories, registers or other such information storage,transmission or display devices. In addition, the term “plurality” maybe used throughout the specification to describe two or more components,devices, elements, parameters and the like.

Embodiments of the invention may manipulate data representations ofreal-world objects and entities such as underground geological features,including faults and other features. The data may be generated bytomographic scanning, as discussed in reference to FIG. 14, e.g.,received by for example a receiver receiving waves generated e.g., by anair gun or explosives, that may be manipulated and stored, e.g., inmemory 150 of FIG. 15, and data such as images representing undergroundfeatures may be presented to a user, e.g., as a visualization on display180 of FIG. 15.

When used herein, a subsurface image or model may refer to acomputer-representation or visualization of actual geological featuressuch as horizons and faults that exist in the real world. Some featureswhen represented in a computing device may be approximations orestimates of a real world feature, or a virtual or idealized feature,such as an idealized horizon as produced in a u_(τ)v_(τ)t_(τ)-transform.A model, or a model representing subsurface features or the location ofthose features, is typically an estimate or a “model”, which mayapproximate or estimate the physical subsurface structure being modeledwith more or less accuracy.

It will thus be seen that the objects set forth above, among those madeapparent from the preceding description, are efficiently attained and,because certain changes may be made in carrying out the above method andin the construction(s) set forth without departing from the spirit andscope of the invention, it is intended that all matter contained in theabove description and shown in the accompanying drawings shall beinterpreted as illustrative and not in a limiting sense.

It is also to be understood that the following claims are intended tocover all of the generic and specific features of the invention hereindescribed and all statements of the scope of the invention which, as amatter of language, might be said to fall therebetween.

1. A system for decompacting a 3D model of the subsurface geology of theEarth at an intermediate restoration time in the past τ, the systemcomprising: one or more processors configured to: receive a 3D model ofpresent-day geometry of the subsurface geology and a measure ofpresent-day porosity experimentally measured within the subsurfacegeology of the Earth, select a value of a restoration time in the past τbefore the present day and after a time an oldest horizon surface in the3D model of the subsurface was deposited, restore the 3D model from thepresent day measured geometry to the predicted past geometry at therestoration time in the past τ using a 3D transformation, and decompactthe vertical dimension of the restored 3D model to elongate verticallengths of geological layers below a horizon layer deposited at therestoration time in the past τ, wherein the vertical lengths areelongated based on a relationship between a depositional porosity of thegeological layers at the time sediment in those layers was deposited,restoration porosity of the geological layers at the restoration time inthe past τ, and the present-day porosity of the geological layersexperimentally measured in the present-day.
 2. The system of claim 1,wherein the one or more processors are configured to decompact thevertical dimension of the restored 3D model by a combination of totaldecompaction corresponding to an increase in porosity from the presentday porosity to the depositional porosity and partial recompactioncorresponding to a partial decrease in the porosity from thedepositional porosity to the restored porosity.
 3. The system of claim1, wherein the restoration porosity is greater than the present-dayporosity and less than the depositional porosity.
 4. The system of claim1, wherein the one or more processors are configured to elongate thevertical length as a measure of height in the vertical dimension.
 5. Thesystem of claim 1, wherein the one or more processors are configured toelongate the vertical length as a measure of geological-time when theparticles of sediment were originally deposited on the Earth's surface.6. The system of claim 1, wherein the one or more processors areconfigured to compute the relationship between the depositionalporosity, the restoration porosity, and the present-day porosity as:$\frac{1 - {\overset{¯}{\varphi}}_{\tau}^{\oplus}}{1 - {\overset{\_}{\varphi}}_{\tau}^{\ominus}},$where ϕ _(τ) ^(⊕)=ψ ₀−ψ and ϕ _(τ) ^(⊖)=ψ ₀−ψ _(p) for all points in therestored 3D model.
 7. The system of claim 1 comprising a probeconfigured to extract material from within the subsurface geology of theEarth or one or more wells and experimentally measure the present-dayporosity of the extracted material.
 8. The system of claim 1, whereinthe one or more processors are configured to use a boundary conditionthat ensures that a top horizon deposited at the restoration time τ is ahorizontal plane in the restored 3D model.
 9. The system of claim 1,wherein the one or more processors are configured to use a boundarycondition that ensures that a direction of change of geological-timewhen the particles of sediment were originally deposited on the Earth'ssurface is vertical in the restored 3D model.
 10. The system of claim 1,wherein the one or more processors are configured to use a boundarycondition that ensures that, for any pair of collocated points onopposite sides of a fault, the two collocated points are decompacted tohave the same coordinate.
 11. The system of claim 1, wherein the one ormore processors are configured to decompact the vertical dimension ofthe restored 3D model by: computing an elongated geological-time in therestored 3D model, transforming the elongated geological-time from therestored 3D model to generate a 1D geological-time in the present-day 3Dmodel, computing 2D paleo-depositional coordinates based on thetransformed geological-time in the present-day 3D model, and performinga 3D transformation comprising the 1D geological-time and 2Dpaleo-depositional coordinates from the present-day 3D model to therestored 3D model that is decompacted based on the elongatedgeological-time.
 12. The system of claim 1, wherein the one or moreprocessors are configured to iteratively decompact the subsurfacelayer-by-layer, starting at the top horizon deposited at the restorationtime τ and ending at the bottom horizon deposited at the depositionaltime.
 13. The system of claim 12, wherein the depositional porosity andthe present-day porosity is independently determined for each geologicallayer of the subsurface.
 14. The system of claim 1, wherein, when thedepositional porosity and the present-day porosity are substantiallyconstant throughout the subsurface geology, the one or more processorsare configured to decompact in one operation over the entire domain ofthe restored 3D model.
 15. A method for decompacting a 3D model of thesubsurface geology of the Earth at an intermediate restoration time inthe past τ, the method comprising: receiving a 3D model of present-daygeometry of the subsurface geology and a measure of present-day porosityexperimentally measured within the subsurface geology of the Earth;selecting a value of a restoration time in the past τ before the presentday and after a time an oldest horizon surface in the 3D model of thesubsurface was deposited; restoring the 3D model from the present daymeasured geometry to the predicted past geometry at the restoration timein the past τ using a 3D transformation; and decompacting the verticaldimension of the restored 3D model to elongate vertical lengths ofgeological layers below a horizon layer deposited at the restorationtime in the past τ, wherein the vertical lengths are elongated based ona relationship between a depositional porosity of the geological layersat the time sediment in those layers was deposited, restoration porosityof the geological layers at the restoration time in the past τ, and thepresent-day porosity of the geological layers experimentally measured inthe present-day.
 16. The method of claim 15 comprising decompacting thevertical dimension of the restored 3D model by a combination of totaldecompaction corresponding to an increase in porosity from the presentday porosity to the depositional porosity and partial recompactioncorresponding to a partial decrease in the porosity from thedepositional porosity to the restored porosity.
 17. The method of claim15, wherein the restoration porosity is greater than the present-dayporosity and less than the depositional porosity.
 18. The method ofclaim 15 comprising elongating the vertical length as a measure ofheight in the vertical dimension.
 19. The method of claim 15 comprisingelongating the vertical length as a measure of geological-time when theparticles of sediment were originally deposited on the Earth's surface.20. The method of claim 15 comprising computing the relationship betweenthe depositional porosity, the restoration porosity, and the present-dayporosity as:$\frac{1 - {\overset{¯}{\varphi}}_{\tau}^{\oplus}}{1 - {\overset{\_}{\varphi}}_{\tau}^{\ominus}},$where ϕ _(τ) ^(⊕)=ψ ₀−ψ and ϕ _(τ) ^(⊖)=ψ ₀−ψ _(p) for all points in therestored 3D model.
 21. The method of claim 15 comprising extractingmaterial from within the subsurface geology of the Earth or one or morewells and experimentally measuring the present-day porosity of theextracted material.
 22. The method of claim 15 comprising applying aboundary condition that ensures that a top horizon deposited at therestoration time τ is a horizontal plane in the restored 3D model. 23.The method of claim 15 comprising applying a boundary condition thatensures that a direction of change of geological-time when the particlesof sediment were originally deposited on the Earth's surface is verticalin the restored 3D model.
 24. The method of claim 15 comprising applyinga boundary condition that ensures that, for any pair of collocatedpoints on opposite sides of a fault, the two collocated points aredecompacted to have the same coordinate.
 25. The method of claim 15comprising decompacting the vertical dimension of the restored 3D modelby: computing an elongated geological-time in the restored 3D model;transforming the elongated geological-time from the restored 3D model togenerate a 1D geological-time in the present-day 3D model; computing 2Dpaleo-depositional coordinates based on the transformed geological-timein the present-day 3D model; and performing a 3D transformationcomprising the 1D geological-time and 2D paleo-depositional coordinatesfrom the present-day 3D model to the restored 3D model that isdecompacted based on the elongated geological-time.
 26. The method ofclaim 15 comprising iteratively decompacting the subsurfacelayer-by-layer, starting at the top horizon deposited at the restorationtime τ and ending at the bottom horizon deposited at the depositionaltime.
 27. The method of claim 26 comprising independently determiningthe depositional porosity and the present-day porosity for eachgeological layer of the subsurface.
 28. The method of claim 15comprising, when the depositional porosity and the present-day porosityare substantially constant throughout the subsurface geology,decompacting in one operation over the entire domain of the restored 3Dmodel.